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Preface 


\ 


In  the  High  Power  Branch  of  the  Air  Force  Aero  Propul- 
sion Laboratory  there  is  a need  for  a reliable  and  light- 
weight pulse-forming  network.  This  need  brought  about  a 
thesis  topic.  In  the  topic  it  was  proposed  that  a student 
develop  a mathematical  model  of  the  network  and  a parameter 
estimation  technique.  The  technique  is  used  to  estimate 
the  inductor  and  capacitor  values  needed  for  the  network  to 
produce  a desired  rectangular  pulse. 

I chose  this  thesis  topic  because  I have  an  under- 
graduate background  in  power  systems  and  a desire  to  better 
understand  these  systems.  In  this  thesis  I not  only  develop 
the  mathematical  model  and  estimation  technique,  but  also 
implement  them  on  the  computer. 

I would  like  to  thank  Dr.  Fredrick  Brockhurst  for  his 
guidance  throughout  all  stages  of  the  thesis  effort  and  Dr. 
J.  Gary  Reid  for  his  assistance  in  developing  the  parameter 
estimation  algorithm.  I would  also  like  to  thank  my  wife 
for  transcribing  the  rough  draft  and  her  tolerance  of  me 
during  the  writing. 


Larry  W . Vannoy 
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Abstract 
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A line-type,  voltage-fed,  pulse-forming  network,  called 
the  "E-line,"  is  studied.  The  energy  to  be  delivered  to  the 
load  is  stored  in  a source  capacitor  at  the  front  of  the  net- 
work. Storing  all  the  energy  in  the  source  capacitor  leads 
to  a more  reliable  and  lightweight  system  than  is  achievable 
by  a conventional  network  where  all  the  capacitors  store  en- 
ergy. The  network  is  to  deliver  100-200  rectangular  pulses 
per  second  to  the  load  where  each  pulse  is  20-30  microsec- 
onds in  duration  and  has  a pulse  height  of  30,000  volts. 

A computer  program  is  developed  which  models  the  dis- 
charge of  the  network  and  is  used  to  analyze  the  output 
pulse  shape  for  different  inductor  and  capacitor  values 
within  the  network.  Another  computer  program  is  developed 
which  estimates  the  inductor  and  capacitor  values  needed  for 
the  network  to  give  a close  approximation  to  the  desired 
rectangular  pulse.  The  simulation  program  models  the  dis- 
charge of  the  actual  network  very  well.  The  estimation 
program  as  applied  in  this  thesis  does  a poor  job  of  pre- 
dicting the  inductor  and  capacitor  values  needed.  The 
desired  rectangular  pulse  is  not  achieved.  Recommendations 
are  included  which  may  improve  the  simulation  program,  esti- 
mation program,  and  the  actual  pulse  shape. 


DIGITAL  SIMULATION  AND  PARAMETER  ESTIMATION 


TECHNIQUES  FOR  THE  E-LINE 
PULSE -FORMING  NETWORK 

I . Introduction 


Background 

In  the  1940 's  one  of  the  technical  problems  facing  the 
scientists  and  engineers  developing  the  microwave  radar  for 
the  armed  services  was  devising  a pulse  generator  capable  of 
delivering  high  power,  short  duration,  and  high  recurrence 
frequency  pulses.  Today,  the  number  of  pulse  generator  ap- 
plications has  grown,  but  their  use  in  radar  systems  stays 
at  the  top  of  the  list.  Along  with  some  of  the  new  applica- 
tions has  come  the  need  for  pulse  generators  which  are  more 
efficient,  lighter  weight,  and  have  higher  power. 

Pulse  generators  depend  on  the  storage  of  electrical 
energy  either  in  an  electrostatic  field  or  in  a magnetic 
field.  Pulse  generators  logically  fall  into  two  major  cat- 
egories, (1)  those  in  which  only  a small  fraction  of  the 
stored  electrical  energy  is  discharged  into  the  load  during 
a pulse,  and  (2)  those  in  which  all  of  the  stored  energy  is 
discharged  during  each  pulse.  These  two  basic  cecegories 
of  pulsers  are  referred  to  as  (1)  "hard-tube  pulsers"  and 
(2)  "line-type  pulsers." 

Line -Type  Pulsers . Line-type  pulsers  are  referred  to 
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as  "line-type"  because  their  energy  storage  device  is  essen- 
tially a lumped-constant  transmission  line.  The  energy  stor- 
age element  serves  not  only  as  a source  of  electrical  energy 
during  the  pulse,  but  also  as  the  pulse  shaping  element, 
thus  it  has  become  commonly  known  as  a "pulse-forming  net- 
work," (PFN).  There  are  essentially  two  classes  of  pulse- 
forming networks,  namely,  those  in  which  the  energy  for  the 
pulse  is  stored  in  an  electrostatic  field  in  the  amount, 

£CV  , called  "voltage-fed  networks,"  and  those  in  which  the 
energy  is  stored  in  a magnetic  field  in  the  amount,  ill/', 
called  "current-fed  networks." 

The  Voltage-Fed  Pulse -Forming  Network.  The  voltage-fed 
pulse-forming  network  in  Figure  1 has  all  its  energy  stored 
in  the  source  capacitor  at  the  front  of  the  network.  This 
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where  the  isolating  element  isolates  the  power  supply 
from  the  current  transients  of  the  discharging  circuit 
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Figure  1.  The  line-type  pulser  on  which 
this  thesis  is  based. 
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where  the  isolating  element  isolates  the  power  supply 
from  the  current  transients  of  the  discharging  circuit 


Figure  2.  The  conventional  line -type  pulser, 


thesis  is  based  on  this  network.  When  the  switch  is  closed, 


the  energy  is  discharged  into  the  load.  Conventional  line- 


type  pulsers  store  energy  in  the  voltage-fed  pulse-forming 


network  by  charging  all  the  capacitive  elements  within  the 
network  (see  Figure  2).  Storing  all  the  energy  in  the 


source  capacitor  leads  to  a much  more  lightweight  and  reli- 


able system  than  is  achievable  by  a conventional  network 


where  all  the  capacitors  store  energy. 


The  reliable,  lightweight  system  has  one  drawback; 


since  the  source  capacitor  must  hold  all  the  network  energy, 


its  value  must  be  much  larger  than  would  be  necessary  if 


the  energy  had  been  shared  among  all  the  capacitive  ele- 


ments of  the  network.  Since  the  values  of  the  capacitances 


and  inductances  of  the  network  determine  the  energy  transfer 


time,  it  is  desirable  to  keep  these  values  small  if  a short 


pulse  width  is  desired.  In  this  thesis  the  lightweight 
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system  of  Figure  1 is  explored  with  an  attempt  to  get  around 
the  pulse  width  problem. 

Problem  Statement 

There  is  a need  for  a mathematical  model  of  the  dis- 
charging circuit  in  Figure  1 which  can  be  implemented  on  the 
computer.  The  model  is  needed  to  analyze  the  output  voltage 
pulse  for  different  sets  of  inductor  and  capacitor  values 
within  the  network.  It  is  also  desirable  to  have  a tech- 
nique for  finding  the  set  of  inductor  and  capacitor  values 
which  will  give  the  desired  output  voltage  pulse. 

In  this  thesis  a computer  program  is  presented  which 
models  the  discharge  of  the  pulse-forming  network  along  with 
a computer  program  which  makes  an  attempt  at  estimating  the 
set  of  inductor  and  capacitor  values  which  will  give  as 
close  as  possible  the  desired  output  voltage  pulse. 

Outline 

Chapter  II  gives  a brief  discussion  on  transmission 
line  theory  and  how  pulse -forming  networks  were  developed 
from  the  theory.  The  particular  network  to  be  explored  in 
this  thesis,  the  "E-line,"  is  also  presented  with  specifica- 
tions to  be  met  by  the  network.  Chapter  III  develops  a 
state  variable  model  of  the  E-line  needed  in  the  network 
simulation  and  the  inductor-capacitor  estimation  programs. 
Chapter  IV  presents  the  simulation  and  estimation  algorithms 
which  are  used  to  form  the  network  simulation  and  parameter 
estimation  computer  programs.  The  programs  can  be  found  in 


Appendix  A and  Appendix  B,  respectively.  In  Chapter  V,  out- 
put voltage  pulses  for  different  sets  of  inductor  and  capac- 
itor values  are  analyzed.  The  final  chapter  contains  con- 
clusions and  recommendations. 

There  are  two  appendices.  Appendix  A contains  a flow- 
chart, operating  instructions,  and  a computer  listing  of  the 
network  simulation  program.  Appendix  B contains  a flow- 
chart, operating  instructions,  and  a computer  listing  of 
the  inductor-capacitor  estimation  program. 
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The  pulse-forming  network,  here  after  referred  to  as 
PFN,  was  developed  from  transmission  line  theory,  in  partic- 
ular from  the  open-ended  lossless  transmission  line.  Al- 
though a PFN  cannot  exactly  simulate  a transmission  line,  a 
close  approximation  to  a transmission  line  is  possible  if 
circuit  elements  of  the  correct  value  are  present. 

This  chapter  begins  by  giving  a brief  discussion  on 
transmission  line  theory.  This  is  followed  by  a section  on 
networks  derived  from  a transmission  line,  and  concludes  by 
discussing  the  PFN  of  particular  interest  in  this  thesis, 
the  "E-line."  Much  of  the  information  in  this  chapter  has 
been  derived  from  the  book,  Pulse  Generators  (Ref  7) . 

Transmission  Line  Theory 

The  transmission  line  that  pulse -forming  networks  try 
to  model  is  the  lossless  open-ended  transmission  line.  The 
nature  of  the  transient  voltage  and  current  produced  by  the 
discharge  of  such  a transmission  line  into  a resistance  load 
(see  Figure  3)  may  be  studied  in  terms  of  its  a-c  impedance 
function. 

The  a-c  impedance  function  of  a lossless  open-ended 
transmission  line  is  given  by  elementary  transmission  line 
theory  as« 

Z = ZQ  coth  j(j)6  (1) 


Figure  3-  Circuit  diagram  of  transmission  line 
discharging  into  a resistive  load 
(Ref  7:177). 


where 


a-c  impedance 

characteristic  impedance  of  the  line 
one  way  transmission  time  of  the  line 


If  = ZQ  , that  is,  if  the  line  is  matched  to  the  load, 
the  current  consists  of  a single  rectangular  pulse  of  ampli- 
tude 1^  = V q/2Z Q and  duration  t = 26  , where  1^  is  the 

transmission  line  current  and  VQ  is  the  voltage  on  the 
transmission  line.  Voltage  and  current  pulses  for  R^  = 

Z0  » Rj  = 2Zq  , and  R^  = (1/2)ZQ  are  shown  in  Figure  4. 

A series  of  steps  is  introduced  into  the  transient  dis- 
charge when  the  line  does  not  match  the  load.  These  steps 


Switch 


Lossless  open-ended 
transmission  line  of 
characteristic  imped 
ance  Z_ . 


where : 

i - length  of  the  transmission  line 
RL  = load  resistor 


Figure  4.  Current  and  voltage  pulses  for  a lossless 
transmission  line  discharging  into  a 
resistive  load.  The  solid  and  broken 
lines  represent  the  voltage  and  current 
pulses,  respectively  (Ref  7sl78). 


are  all  positive  when  is  larger  than  ZQ  and  alternate  in 
sign  when  R^  is  smaller  than  ZQ.  The  steps  are  due  to  re- 
flections caused  at  the  terminals  of  the  line  by  mismatching 
the  load  resistance.  "These  reflections  traverse  the  line 
to  the  open  end  in  time  6 , and  are  completely  reflected 
there,  and  travel  back  to  the  load  end  in  a total  time  26, 
where  they  appear  as  positive  or  negative  steps  depending 
upon  the  mismatched  ratio"  (Ref  7tl78-l79).  The  reflections 
continue,  with  diminishing  amplitude,  until  all  the  energy 
initially  stored  in  the  line  is  dissipated  in  the  resistive 
load. 

One  might  ask,  "Why  not  use  a transmission  line  rather 
than  a PFN  if  a rectangular  pulse  is  desired?"  A 
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transmission  line  is  not  used  for  producing  rectangular 
pulses  because  the  diameter  of  the  transmission  line  needed 
for  high  voltage  applications  and  the  length  the  line  would 
have  to  be  for  the  correct  pulse  duration  is  not  practical 
for  most  applications. 

Networks  Derived  From  a Transmission  Line 

The  development  of  PFN's  which  simulate  a transmission 
line  is  a mathematical  problem  in  network  synthesis.  At 
present  there  are  several  PFN's  that  have  been  developed 
(see  Figure  5)*  It  should  be  noted  that,  although  a PFN  can 
simulate  a transmission  line,  no  network  having  a finite 
number  of  lumped  parameter  elements  can  exactly  simulate  a 
transmission  line  which  in  reality  has  distributed  parame- 
ters. As  the  number  of  lumped  parameter  elements  in  a given 
PFN  increases,  the  degree  of  approximation  will  improve. 

Pulse  forming  networks  are  essentially  ladder  networks. 
The  sections  or  windows  of  the  network  are  called  meshes. 

The  voltage  drop  across  and  current  in  any  mesh  are  not  only 
dependent  on  the  elements  within  the  mesh  but  also  the  volt- 
age and  current  in  all  the  other  meshes  taken  separately  and 
together  as  a group.  The  PFN  will  not  produce  a rectangular 
pulse  if  there  are  not  enough  meshes  in  the  network  or  the 
elements  within  each  mesh  have  the  wrong  value;  for  example, 
the  network  pulse  may  rise  too  slowly  or  have  excessive 
overshoots  at  the  beginning  of  the  pulse,  may  ripple  through- 
out the  pulse  width,  and  may  fall  too  slowly  or  have  over- 
shoots at  the  end  of  the  pulse.  It  is  part  of  this  thesis 
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to  develop  a technique  to  find  an  approximate  value  for  the 
elements  within  each  mesh  of  a particular  PFN  as  the  number 
of  meshes  changes. 

The  E-Line  Pulse -Forming  Network 

The  rest  of  this  thesis  will  be  based  on  the  E-line 
pulse-forming  network  (see  Figure  6).  As  mentioned  in  the 
introduction,  this  network  differs  from  the  type  used  by  a 
conventional  line-type  pulser  scheme  in  that  the  energy  to 
be  delivered  to  the  load  has  all  been  stored  in  the  first 
capacitor  of  the  network  rather  than  charging  up  all  the 
capacitors  of  the  network.  Storing  all  the  energy  in  the 
first  capacitor,  rather  than  storing  some  of  the  energy  in 
all  of  the  capacitors  of  the  network,  causes  the  first  ca- 
pacitor of  the  network  to  be  several  times  larger  than  would 
be  required  if  the  energy  storage  was  shared.  The  E-line 
of  interest  has  the  following  specif icationsj 

1.  It  produces  100-200  rectangular  pulses  per  second. 

2.  Each  pulse  has  a 20-30  microsecond  duration. 

3.  All  possible  energy  is  transferred  from  the  source 
capacitor  to  the  load. 

U.  The  source  capacitor  has  been  charged  to  30,000 
volts  and  contains  60,000  joules  of  energy 
initially. 

5.  Voltage  pulse  height  is  on  the  order  of  30,000 
volts . 

6.  There  are  to  be  few  oscillations  after  the  pulse. 

7.  The  charging  capacitor  is  charged  between  pulses. 

8.  The  load  is  purely  resistive. 

9.  At  initial  start  up,  initial  conditions  are  zero 
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on  all  circuit  elements  except  the  charging  ca- 
pacitor. 


10.  The  E-line  will  contain  5~7  meshes. 

Transmitting  most  of  the  energy  from  the  source  capacitor 
to  the  load  is  the  most  important  specification  in  the  list. 

In  any  pulse -forming  network  it  is  desirable  to  have 
the  characteristic  impedance  match  the  load  impedance  so  no 
reflections  occur  in  the  line.  However,  this  is  not  always 
possible.  The  characteristic  impedance  of  the  E-line  ana- 
lyzed in  this  thesis  in  fact  does  not  match  the  load  imped- 
ance and,  therefore,  reflections  add  another  dimension  to 
the  problem  of  getting  a rectangular  pulse.  Hopefully,  the 
E-line  element  approximation  routine  to  be  discussed  later 
in  the  text  will  minimize  the  reflections  within  the  line. 

Basic  transmission  line  theory  has  been  presented  in 
this  chapter  along  with  some  pulse  forming  networks  which 
approximate  the  open-ended  lossless  transmission  line.  It 
should  be  remembered  that  a network  composed  of  a finite 
number  of  lumped  parameter  elements  cannot  exactly  simulate 
the  transmission  line  which  in  reality  has  distributed 
parameters . 
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III.  A State  Variable  Model  of  the  E-Line 

State  equations,  a definition  of  state,  and  the  solu- 
tion to  the  state  equations  are  presented  in  this  chapter 
along  with  the  set  of  state  equations  that  describes  the  E- 
line  PFN.  In  the  next  chapter  the  state  equations  will  be 
used  to  build  an  analog  model  of  the  E-line  and  a parameter 
estimation  algorithm  which  estimates  the  correct  capacitor 
and  inductor  values  to  give  a rectangular  pulse. 

State  and  State  Equations 

"The  state  of  a system  is  a mathematical  structure 

containing  a set  of  n variables  X^( t) , X2( t) , ' * * , X^(  t) , * ’ ' , 

X ft),  called  state  variables,  such  that  the  initial  values 
n 

Xi(to)  of  this  set  and  the  system  inputs  U^(t)  are  suffi- 
cient to  uniquely  describe  the  system's  future  response  for 
t > t " (Ref  1:23-24).  The  E-line  of  interest  has  no  in- 
puts, all  initial  conditions  are  zero  except  for  the  initial 
voltage  on  the  source  capacitor,  and  physical  variables  are 
chosen  as  the  state  variables. 

The  mathematical  structure  mentioned  in  the  above  quo- 
tation (Ref  1:23-24)  is  a set  of  first  order  differential 
equations  called  state  equations.  These  equations,  with  the 
initial  condition  on  the  source  capacitor,  are  able  to  com- 
pletely describe  the  PFN’s  output  time  response. 

Physical  Variables . The  selection  of  physical  vari- 
ables as  the  state  variables  is  based  upon  the  energy  stor- 
age elements  in  the  PFN  (see  Table  I).  The  variables  must 
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TABLE  I 

Energy  Storage  Elements  (Ref  1:25) 


Element 

Energy 

Physical  variable 

Capacitor  C 

Cv2 

2 

Voltage  v 

Inductor  L 

Li2 

2 

Current  i 

be  independent,  that  is,  they  cannot  be  expressed  in  terms 
of  the  remaining  physical  variables.  Since  there  is  a 
source  capacitor  at  the  front  of  the  E-line  and  two  storage 
elements  in  each  mesh,  then  there  are  2n  + 1 physical  vari- 
ables where  n is  the  number  of  meshes. 

State  Equations . Since  there  are  2n  + 1 state  vari- 
ables there  must  be  2n  + 1 state  equations.  Since  this  is 
a first  attempt  at  deriving  state  equations  for  the  E-line, 
mutual  inductance  in  the  network  has  been  neglected.  For 
simplification,  the  derivative  of  a variable  is  signified 
by  putting  a dot  over  the  variable . 

The  state  equations  for  the  E-line  are  the  loop-current 
and  node-voltage  equations  of  each  mesh  (see  Figure  7).  For 
example,  the  voltage  across  CS  is 


.rt  ( 
CS  J 


ILldt  + vcs(°> 


or  the  state  equation  is 


v = t 
VCS  CS  iLl 
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Notice  the  initial  condition,  Vcs(0). 

The  current  in  inductor  is 

IL1  = J ^VCS  " R1IL1  ‘ VCl^dt 

o 

or  the  state  equation  is 


The  state  equations  for  the  rest  of  the  L's  and  C’s  then 


become 


where  n is  the  number  of  meshes 


The  state  equations  can  be  expressed  in  matrix  form  as 


where 


an  mxl  column  vector  of 
state  variable  derivatives 


an  mxm  plant 
coefficient  matrix 


an  mxl  column  vector 
of  state  variables 


m = 2n  + 1 

n = number  of  meshes 

The  output  equation  can  then  be  expressed  as 


Y = CTX 


where 


in 

C = [C.  C2  . . . C l an  lxm  output 
1 row  vector 


The  E-line's  state  equations  expressed  in  matrix  form 
becomes  Equation  (18)  shown  at  the  top  of  the  next  page. 
Since  the  voltage  on  the  load  is  the  primary  interest,  the 
output  equation  becomes 


Y = [000 


(17) 


which  is  a scalar. 


i The  state  variables  and  state  equations  developed  here 

are  not  the  only  ones  that  can  describe  the  E-line j however, 
these  variables  and  equations  are  more  readily  understandable 
and  logically  derived  than  some  of  the  other  forms . If 
there  had  been  any  inputs  on  the  PFN  ( inputs  do  not  include 
initial  conditions),  the  state  and  matrix  form  of  the  state 
equations  would  have  been  more  complicated}  for  information 
on  this  situation  see  Ref  It  14-28. 


Solving  the  Output  Equation . For  a system  which  has 
only  initial  conditions  and  the  initial  time,  tQ,  is  taken 
as  zero,  the  solution  to  the  output  equation  is 
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u 


o 


Y(t)  = CT  exp  [At]X(0) 


(19) 


where 


C1  = the  lxm  output  row  vector 
A = the  mxm  coefficient  matrix 

exp  [At]  , I t Iff  + ♦.! AU2  + . . . ifiil 

= an  mxm  matrix  called  the  "state  transition 
matrix" 

I = the  identity  matrix 


X(0)  = 


X^o) 

x2(0) 

x3(0) 


vo) 


an  mxl  vector  containing 
the  initial  conditions  of 
the  state  variables 


When  this  equation  is  applied  to  the  E-line,  the  CT  vector 
is  given  in  Eq  (18),  the  A matrix  is  given  in  Eq  (19)»  and 
the  initial  condition  vector,  X(0),  is 


X(0)  = 


Vcs(°) 


0 

p 

Note  that  the  only  initial  condition  in  the  network  is  the 
initial  voltage  on  the  source  capacitor,  and  the  solution  to 
the  output  equation  is  actually  the  voltage  across  the  load. 
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The  next  two  chapters  will  use  this  solution  in  an  attempt 
to  find  the  inductor  and  capacitor  values  needed  for  the 
E-line  to  produce  a rectangular  voltage  pulse  across  the 
load . 

There  are  many  methods,  other  than  the  one  above,  for 
evaluating  the  state  transition  matrix.  If  the  A matrix 
becomes  larger  than  3x3,  evaluating  the  state  transition 
matrix  by  hand  can  be  quite  a task;  therefore,  many  com- 
puter routines  have  been  developed.  The  subroutine  MEXP 
(Ref  6; 1.13),  presented  in  the  next  chapter,  is  used  in  the 
inductor  and  capacitor  approximation  algorithm  to  solve  for 
the  state  transition  matrix. 

MEXP  (Ref  6:1.13)  is  a subroutine  developed  at  the 
University  of  Connecticut's  School  of  Engineering.  The 
subroutine  solves  the  state  transition  matrix  using  a 
Chebyshev  approximation  rather  than  the  series  expansion 
technique  presented  here . 

The  information  in  this  chapter  only  touches  the  sur- 
face of  the  state  variable  and  state  equation  concepts. 

Most  of  the  material  in  this. chapter  came  from  the  book, 
Linear  Control  System  Analysis  and  Design  (Ref  1:14-85), 
however,  this  material  can  be  found  in  almost  any  intro- 
ductory book  on  modem  control  theory.  This  material  will 
be  used  extensively  in  the  next  chapter. 

o 
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IV . Digital  Simulation  and  Parameter  Estimation 
Algorithms  for  the  E-Line  Pulse-Forming  Network 


( 


The  first  section  in  this  chapter  begins  by  building 
an  analog  model  of  the  E-line.  Then  a simulation  algorithm 
will  be  developed  using  a computer  subroutine  which  approx- 
imates an  analog  computer,  RKFOUR  (see  Appendix  A).  The 
computer  program  formed  from  the  simulation  algorithm  will 
be  used  in  the  next  chapter  to  analyze  the  PFN  as  the  values 
of  the  inductors  and  capacitors  within  the  network  change. 

In  the  second  section  a parameter  estimator  algorithm  is 
developed  which  estimates  the  inductor  and  capacitor  values 
needed  in  the  E-line  for  the  network  to  produce  a rectangu- 
lar output  voltage  pulse.  The  estimation  algorithm  uses  a 
function  minimization  subroutine,  ZXMIN,  from  the  Interna- 
tional Mathematics  and  Statistical  Libraries,  (Ref  5:ZXMIN1- 
ZXMIN3)  . The  computer  programs  formed  from  the  simulation 
and  parameter  estimation  algorithms  can  be  found  in  Appendix 
A and  Appendix  B,  respectively. 
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Analog-Digital  Simulation  Algorithm  of  the  E-Line 

As  in  the  state  equations,  loop-current  and  node -volt- 
age equations  will  be  used  to  develop  an  analog  model  of 
the  PFN.  However,  first  order  integral  equations  rather 
than  first  order  differential  equations  are  needed. 

Analog  integrator  and  summing  elements  are  shown  in 
Figure  8.  When  making  an  analog  model  of  the  E-line,  the 
input  to  a particular  integrator  is  either  the  voltage 
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Integrator 

Z< 


where : Z^input,  G=gain,  Z2=output 


Summer 


where:  ,U2,U^= inputs,  U^=output 


Figure  8.  Analog  Integrating  and  Summing  Elements 


across  an  inductor  or  the  current  through  a capacitor.  The 
output  of  a particular  integrator  is  either  the  current 
through  an  inductor  or  the  voltage  across  a capacitor. 

In  Figure  7 the  voltage  across  the  source  capacitor  is 
given  by 


(20) 


The  current  in  the  inductor,  L^t  is 


■k.f 


(21) 


When  the  equations  are  put  in  analog  diagram  form,  the 
diagram  turns  out  to  be  that  of  Figure  9. 

The  subroutine  RKFOUR  (see  Appendix  A)  is  a computer 
routine  that  simulates  an  analog  computer.  The  equations 
used  to  form  the  analog  diagram  of  Figure  9 will  be  put  into 
the  program  and  the  performance  of  the  E-line  will  be  ana- 
lyzed. B.D.  Weathers  developed  RKFOUR  while  at  the  Univer- 
sity of  Missouri,  Columbia  Campus. 

The  program  that  calls  RKFOUR  has  been  set  up  to  handle 
varying  pulse  rates,  varying  numbers  of  meshes,  and  varying 
initial  voltages  on  the  source  capacitor  so  that  E-line 
specifications  (1),  (4),  and  (9)  on  page  11  can  be  met. 
However,  getting  the  energy  transferred  to  the  load,  having 
a rectangular  pulse  of  correct  duration,  and  having  a 30,000 
volt  pulse  height  are  dependent  on  choosing  the  correct  L's 
and  C's  in  the  network. 

A switching  technique  is  used  to  help  get  the  energy 
from  the  source  capacitor  to  the  load.  Some  of  the  energy 
in  the  network  starts  going  back  into  the  source  capacitor, 
CS,  if  while  CS  is  discharging  its  energy  into  the  network, 
the  current  in  inductor  goes  negative . The  recharging 
of  CS  is  due  to  reflections  in  the  network.  The  reflections 
cause  the  discharge  of  CS  to  be  slow  and  a 20-30  microsecond 
pulse  width  is  not  achievable.  To  at  least  prevent  the  re- 
flections from  getting  back  to  CS,  the  switch  between  it  and 
opens  when  the  current  in  inductor  starts  to  go  nega- 
tive. When  the  voltage  across  CS  becomes  larger  than  the 


Figure  9.  An  Analog  Model  of  the  E-Line 


voltage  across  capacitor  , and  the  current  would  normally 
start  going  positive  again,  the  switch  is  closed.  This 
switching  shortens  the  amount  of  time  needed  to  discharge 
CS  into  the  load.  When  the  voltage  across  CS  reaches  a pre- 
determined value  and  the  current  again  starts  to  go  negative 
the  capacitor  is  removed  and  recharged  for  the  next  pulse 
while  the  rest  of  the  network  is  allowed  to  run.  The  switch 
ing  occurs  after  each  pulse  during  the  entire  pulse  train. 
The  switching  has  also  been  programmed  into  the  main  program 
that  calls  RKFOUR . 

Appendix  A contains  a flowchart,  the  operating  instruc- 
tions, and  a computer  listing  of  RKFOUR.  In  the  operating 
instructions  and  computer  printout  there  should  be  enough 
information  to  understand  how  RKFOUR  works  and  how  it  is 
used . 

Parameter  Estimation  Algorithm 

It  would  be  nice  if  a computer  routine,  rather  than 
trial  and  error,  could  be  used  to  find  the  inductor  and 
capacitor  values  that  will,  as  close  as  possible,  approxi- 
mate the  desired  rectangular  voltage  pulse  at  the  output. 

In  this  section  a computer  program  is  developed  that  uses 
a least  squares  fit  objective  function  in  a function  mini- 
mizing subroutine,  ZXMIN  (Ref  51ZXMINI-ZXMIN3) » to  determine 
the  L's  and  C’s  needed  to  give  the  rectangular  pulse.  Since 
initial  methodology  is  being  developed  here,  the  estimation 
algorithm  only  tries  to  approximate  a rectangular  voltage 
pulse  for  the  first  pulse  in  a series  of  pulses.  The  first 


pulse  is  approximated  because  the  simplest  case  is  when  all 
initial  conditions  are  zero  except  for  Vcs(0).  Hopefully, 
if  a rectangular  pulse  is  achieved  for  the  first  pulse  then 
the  rest  of  the  pulses  in  the  series  will  also  be  rectangu- 
lar. A more  desirable  pulse  to  approximate  would  be  a pulse 
later  in  the  pulse  train  after  the  prepulse  initial  condi- 


tions have  reached  steady  state 


ZXMIN  is  a subroutine  which  minimizes  a function  of  n 


variables  using  a quasi-Newton  algorithm.  The  function  to 


be  minimized  in  this  application  is  the  least  squares  fit 
objective  function  and  the  variables  are  the  inductor  and 
capacitor  values  in  the  network,  ZXMIN  calls  the  subroutine 
that  contains  the  least  squares  fit  objective  function, 
FUNCT,  and  by  trying  different  L's  and  C's  minimizes  the 
function.  A printout  of  the  operating  instructions  for 


ZXMIN  is  in  Appendix  B 


The  E-line’s  desired  output  voltage  pulse  is  shown  in 
Figure  10  with  what  the  actual  output  voltage  pulse  might 
look  like.  The  pulse  does  not  start  at  time  TS  = 0 be- 
cause there  is  a finite  time  before  the  voltage  on  the 
source  capacitor  can  get  to  the  load.  It  is  desirable  for 
the  actual  pulse  to  not  only  approximate  the  desired  pulse 
but  go  to  zero  after  time  Tp. 

The  least  squares  fit  objective  function  is  defined  as 


Figure  10.  Desired  and  Actual  Output  Pulses 


where 

= weighting  factor 

When  F is  a minimum  the  calculated  pulse  approximates  the 
desired  pulse  as  closely  as  possible  while  keeping  the  os- 
cillations after  the  pulse  to  a minimum.  The  oscillations 
after  the  pulse  are  due  to  continued  reflections  within  the 
E-line.  The  weighting  factor,  K,  determines  how  desirable 
it  is  to  keep  the  postpulse  oscillations  small.  If  K1  is 
small,  0 < <0.5*  then  approximating  the  pulse  is  more 

desirable  than  the  postpulse  oscillations.  If  is  large, 
0.5  < < 1.0,  then  the  postpulse  oscillations  are  at  least 

as  desirable  as  approximating  the  pulse. 

Y(t)  in  Eq  (29)  is  the  output  voltage,  of  Figure  ?. 


* 


' 

The  object  of  minimizing  Eq  (29)  is  to  make  the  voltage  ap- 
proximate the  desired  rectangular  pulse.  Recall  from  Eq  (19) 
of  Chapter  III  that  Y(t)  is  also  known  as  the  solution  to 
the  output  equation,  i.e., 

Y( t)  - VRL  = CT  exp  [At]X(0)  (30) 

where,  for  the  E-line: 

C = the  ixm  output  row  vector  [00* ••01] 
exp  [At]  = the  mxm  state  transition  matrix 
X(0)  = the  mxl  initial  condition  vector 

= VCS^0^" 

0 


0 


Therefore,  the  solution  to  the  output  equation  can  be,  and 
is,  used  for  Y(t)  in  the  least  squares  fit  objective  func- 
tion. 


The  integrals  in  Eq  (29)  are  done  numerically  using 
the  trapezoidal  rule  (Ref  4:144-148).  According  to  the 
trapezoidal  rule,  to  evaluate  the  integral 


I 


f(  t)dt 


(3D 


one  can  divide  the  interval,  a < t < b,  into  n equal  sub in- 
tervals each  of  width  At  (see  Figure  11),  where 
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Figure  11.  Approximating  an  Integral 


and  the  integral  is  then  approximated  by 


/ f(t)dt  « |i(f0+2fi+2f2+---  + 2fn_2+2fn_i+fn) 


where 

.fo. • • *fL  , are  the  values  of  the  function  at 
0 l 2 n-i  n 

the  points  tQ, t1 ,t2* • *tn_lttn 


As  At  gets  smaller,  the  approximation  to  the  true  integral 
becomes  more  accurate. 

When  applying  the  trapezoidal  rule  to  the  least  squares 


(32) 
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At  = 
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fit  objective  function,  Eq  (29),  the  result  is 


* - Is  < 


TP 


(PH  - Y(TS) ) ^ + 


(r  At 

2 I 

k , TS 


(PH  - Y(i-At)) 


• TS  < 
1=At  1 


') 


3TP 

('  At 

■i  I 

> . TP 


') 


(Y(i-At))2  ] + (Y(3TP))2 


(34) 


Notice  that  J|,  and  ^ must  come  out  to  be  integers. 

When  substituting  Eq  (30)  for  Y(TS),  Y(i*At),  and  Y(3TP) 


Eq  (34)  becomes 


v _ At 
F " 2 


(PH  - CT  exp  [A*TS]X(0))2 


V TP 
At 


(’  At 

• l 

k . _TS 
1 At 


(PH  - CT  exp  [A*i*At]X(0)) 


!) 


^ £ (CT  exp  [A*  i*At]X(0))^ 


+ (CT  exp  [A-3TP]X(0)) 


(35) 


If  the  state  transition  matrix  is  found  for  time  At,  i.e., 


exp  [A.t]  06) 


,k/  Aa.Nk 
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an  mxm  matrix 


the  values  for  exp  [ATS^,  exp  [Ai*At],  and  exp  [A3TP]  can  be 
found  easily  by  bringing  the  value  of  exp  [AAt]  to  the  cor- 
rect power,  i .e . , 


exp  [A’TS] 
exp  [Ai*At] 


exp  [(AAt)  • i] 


(37) 

(38) 


and 


exp  [A*3TP] 


(39) 


The  subroutine  that  contains  the  least  squares  fit  objective 
function,  FUNCT,  calls  the  subroutine  MEXP  (see  page  22) 
which  calculates  the  state  transition  matrix  exp  [A*At]. 

Earlier  in  this  section  it  was  stated  that  the  subrou- 
tine ZXMIN  tries  to  minimize  the  subroutine  that  contains  the 
least  squares  fit  objective  function,  actually  ZXMIN  tries  to 
minimize  Eq  (35)»  an  approximation  to  the  least  squares  fit 
objective  function.  Appendix  B contains  a flow  chart  of 
FUNCT,  operating  instructions  for  ZXMIN  and  MEXP,  a computer 
printout  of  the  main  program  which  calls  ZXMIN,  and  a computer 
printout  of  subroutine  FUNCT. 

In  this  chapter  the  algorithms  needed  to  form  the  ana- 
log approximation  routine  and  the  parameter  estimation  rou- 
tine were  developed.  The  next  chapter  uses  the  programs  to 
compare  sets  of  L’s  and  C’s  and  predict  trends  for  the  L's 
and  C's  to  get  the  desired  rectangular  output  voltage  pulse. 
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V . Analysis  of  the  Simulation ' s Results 

The  analog  approximation  routine  described  in  the 
previous  chapter  has  been  set  up  to  simulate  the  voltage  and 
current  time  responses  for  the  E-line.  The  values  for  the 
inductors  and  capacitors  predicted  by  the  parameter  estima- 
tor are  put  into  the  simulation  program  with  several  other 
sets  of  inductor  and  capacitor  values.  The  output  voltage 
calculated  from  the  simulation  program  is  plotted  on  graphs 
for  the  different  sets  of  inductor  and  capacitor  values. 
Hopefully,  the  graphs  will  aid  in  predicting  a trend  for 
the  inductors  and  capacitors  to  follow  in  achieving  the 
desired  rectangular  pulse. 

The  E-line  for  which  all  data  is  based  has 

1.  5 meshes. 

-4 

2.  Source  capacitance,  CS,  of  1.33x10  farads. 

3.  A pulse  repetition  rate  of  100  pulses  per  second. 

4.  A source  capacitor  initial  voltage  of  30,000 
volts . 

The  desired  pulse  is  delayed  by  5 microseconds,  has  a pulse 
width  of  20  microseconds,  and  has  a pulse  height  of  30,000 
volts,  i.e., 

TS  = 5x10"^  sec 
TP  = 25x10" 6 sec 
PH  = 30,000  volts 

Data  on  several  more  E-line  and  desired  pulse  cases  would 
have  been  collected  but  time  was  not  permitting. 
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L and  C Values  Estimated  by  ZXMIN 

When  the  preceding  values  for  TS,  TP,  and  PH  are  used 
in  the  least  squares  fit  objective  function  and  the  weight- 
ing factor,  K1 , is  chosen  as  0.5»  the  parameter  estimation 
routine  predicts  the  capacitances  and  inductances  given  in 
row  1 of  Table  II.  When  these  parameter  values  are  put  into 
the  E-line,  the  resulting  output  voltage,  for  the  weighting 
factor  chosen,  should  approximate  as  close  as  possible  the 
desired  rectangular  pulse. 

Base  Values  for  L and  C 

In  the  book  Pulse  Generators  (Ref  7*201)  a set  of  in- 
ductor and  capacitor  values  has  been  calculated  which  gives 
an  approximately  rectangular  output  voltage  pulse  if  a con- 
ventional line-type  pulser  charging  scheme  is  used  to  put 
energy  into  the  network  (see  page  3 of  this  thesis).  Al- 
though the  E-line  of  interest  uses  a different  type  charging 
scheme,  the  values  calculated  in  Ref  7 are  used  as  base 
values  for  the  inductors  and  capacitors  in  the  network  (see 
Figure  12) . Although  these  values  are  not  the  exact  values 
needed  by  the  E-line  to  produce  the  desired  rectangular 
pulse,  they  are  a starting  point.  The  base  values  along 
with  multiples  of  the  base  values  are  found  in  Table  II. 
Because  the  last  capacitor  in  the  network,  C^,  is  important 
in  determining  the  time  constant  at  the  load,  should  not 
be  changed;  thus,  stays  the  base  value  in  all  multiples 
of  the  base  values. 
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Figure  12.  Multiply  the  values  of  the  inductances  by 

the  product  RL  (Pulse  Width)  and  the  values 
of  the  capacitances  by  Pulse  Width/RL. 

The  inductances  are  in  henrys  and  the  ca- 
pacitances in  farads.  (Ref  7:201) 


Analysis  of  Results 

The  simulation  program  has  calculated  the  output  volt- 
age for  each  of  the  test  cases  in  Table  II.  Graphs  of  the 
output  voltage  as  a function  of  time  can  be  found  in  Figures 
13,  14,  ani  15-  The  desired  pulse  shape  has  been  drawn  on 
all  the  graphs,  and  the  graphs  have  been  normalized  to  the 
30,000  volt  desired  pulse  height. 

In  Figure  13  it  can  be  seen  that  although  the  inductor 
and  capacitor  values  predicted  by  the  parameter  estimator  do 
not  produce  a curve  that  approximates  the  desired  pulse  very 
well,  there  is  a noticeable  improvement  over  the  curve  re- 
sulting from  the  base  values  of  L and  C.  When  vising  the 
estimator's  values  the  rise  time  is  very  good  and  the  first 
two  peaks  fall  within  the  desired  pulse  width;  however,  the 
fall  time  is  far  too  long  and  the  pulse  height  is  not  nearly 
high  enough.  The  base  values'  voltage  curve  has  a slightly 
higher  first  peak,  but  the  second  peak  falls  outside  the 
desired  pulse  width  and  is  too  short;  also,  the  fall  time 
is  too  long. 

In  Figure  14  the  curves  shift  to  the  right  and  the 
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The  E-line's  output  voltage  as  calculated  by  the  simulation  program  for 
different  sets  of  inductor  and  capacitor  values. 


Figure  14.  The  E-line's  output  voltage  as  calculated  by  the  simulation  program  for 
different  sets  of  inductor  and  capacitor  values. 


Figure  15.  The  E-line’s  output  voltage  as  calculated  by  the  simulation  program  for 
different  sets  of  inductor  and  capacitor  values. 


height  of  the  second  peaks  decrease  as  the  capacitor  values 
in  the  network  increase.  When  the  base  values  for  the  ca- 
pacitor are  cut  in  half  the  corresponding  curve  has  its 
first  two  peaks  at  about  the  same  height  and  within  the 
desired  pulse  width,  but  there  are  three  undesirable  oscil- 
lations outside  the  pulse  and  the  fall  time  is  again  too 
long.  As  the  capacitor  values  are  increased  above  their 
base  value,  their  corresponding  curves  are  definitely  unde- 
sirable in  that  there  is  only  one  major  peak  and  the  fall 
time  is  even  longer. 

In  Figure  15,  as  the  inductor  values  increase,  the 
curves  shift  to  the  right  and  the  height  of  the  first  peak 
decreases.  When  the  base  values  for  the  inductors  are  cut 
in  half,  the  corresponding  curve  has  a very  good  rise  time, 
has  high  first  and  second  peaks,  and  the  fall  time  between 
the  second  and  third  peaks  looks  good.  However,  the  curve's 
overall  fall  time  is  too  long  and  there  are  two  undesirable 
oscillations  outside  the  desired  pulse  width.  As  the  in- 
ductor values  are  increased  it  takes  too  long  for  the  pulse 
to  rise,  there  are  pronounced  oscillations  after  the  pulse, 
and  the  pulse  width  has  to  be  increased  for  most  of  the  en- 
ergy in  the  network  to  get  to  the  load. 

None  of  the  curves  on  Figures  13,  14,  or  15  look  much 
like  the  desired  output  pulse.  When  the  weighting  factor 
is  chosen  as  0.5  it  is  apparent  that  the  parameter  estimator 
did  a poor  job  of  predicting  the  inductances  and  capacitances 
needed  for  the  network  to  approximate  the  desired  pulse. 
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Changing  the  weighting  factor  or  making  the  weighting  factor 
time  varying  may  improve  the  parameter  estimator’s  predic- 
tion. 

Trends 

Although  the  curves  are  somewhat  different  from  the 
desired  pulse,  the  curves  do  display  some  trends.  If  the 
inductor  values  are  decreased  from  their  base  values,  the 
rise  time  and  fall  time  can  be  improved.  If  the  capacitor 
values  are  decreased  from  their  base  value,  except  for  the 
last  capacitor,  the  first  and  second  peaks  can  be  made  to 
have  the  same  height.  Also,  if  the  inductor  values  are  de- 
creased the  first  and  second  peaks  can  both  be  higher.  How- 
ever, if  the  inductor  and  capacitor  values  are  decreased  too 
much  the  output  voltage  will  look  as  though  the  source  ca- 
pacitor is  discharging  directly  into  the  load. 

If  the  trends  mentioned  above  are  incorporated  into  a 
trial  and  error  technique  it  may  be  possible  to  find  the 
inductor  and  capacitor  values  needed.  Trial  and  error 
techniques,  however,  can  be  quite  time  consuming  and  the 
time  might  be  better  spent  on  attempting  to  improve  the 
existing  estimation  program  or  develop  a new  program. 


VI.  Conclusions  and  Recommendations 


Conclusions 

The  network  simulation  program  models  the  discharge, 
switching  action,  and  pulsing  action  of  the  actual  network 
very  well  and  is  indispensable  in  analyzing  the  output  volt- 
age for  different  sets  of  inductor  and  capacitor  values.  The 
program  is  flexible  enough  to  model  a variety  of  changes 
that  may  be  desired  in  the  actual  network. 

The  parameter  estimation  program,  as  applied  in  this 
thesis,  is  not  working  as  well  as  desired.  It  may  not  be 
working  well  because  of  a simple  FORTRAN  error,  but  it  is 
more  likely  due  to  an  error  in  the  estimation  algorithm, 
not  describing  the  weighting  factor  properly,  and/or  prob- 
lems with  the  subroutines  used  by  the  program.  Unless  the 
estimation  program  can  be  improved  or  another  estimation 
program  can  be  developed,  finding  the  inductor  and  capacitor 
values  needed  for  the  E-line  to  approximate  the  desired  out- 
put voltage  pulse  may  involve  considerable  trial  and  error. 
Recommendations  are  made  in  the  next  section  on  what  might 
be  done  to  improve  the  parameter  estimator  or  develop  a 
new  parameter  estimator. 

The  network’s  desired  output  voltage  pulse  was  not 
achieved  with  any  of  the  test  cases  listed  in  Table  II,  but 
the  curves  corresponding  to  the  test  cases,  Figures  13,  14, 
and  15,  displayed  some  trends.  Using  these  trends  as  a 
starting  point,  the  desired  rectangular  pulse  may  be 
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achievable  through  trial  and  error.  However,  it  may  be  the 
case  that  it  is  not  possible  to  approximate  the  desired  rec- 
tangular pulse  very  well  with  the  particular  E-line  pulse - 
forming  network  analyzed.  In  the  following  section,  changes 
that  might  be  made  in  the  network  to  get  the  calculated  out- 
put voltage  pulse  are  recommended. 

Recommendations 

Although  the  simulation  program  seems  to  describe  the 
actual  output  voltage  well,  including  mutual  inductance  in 
the  program  will  improve  its  accuracy.  Mutual  inductance 
has  a significant  impact  on  the  output  voltage  wave  shape 
(Ref  2) . Therefore,  if  the  mutual  inductance  is  included  in 
the  simulation  program,  it  may  be  found  that  the  desired 
rectangular  pulse  has  been  achieved.  If  this  thesis  is 
chosen  as  a starting  point  for  another  thesis,  mutual  in- 
ductance should  definitely  be  included  in  the  simulation 
program . 

As  is,  the  parameter  estimator  program  is  doing  a poor 


job  of  predicting  the  inductor  and  capacitor  values  needed. 
To  improve  the  estimator's  predicting  ability,  one  or  more 
of  the  following  might  be  tried: 


1. 

Include  more  of  the  system's  dynamics,  such  as 
mutual  inductance  between  meshes. 

2. 

Have  the  estimator  try  to  approximate  a whole 
series  of  pulses  in  the  pulse  train. 

3- 

Change  the  value  of  the  weighting  factor  or  make 
the  weighting  factor  time  varying. 

1 4. 

Make  it  possible  for  the  estimation  program  to 
only  change  the  inductor  and  capacitor  values, 

not 

but 
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also  the  resistor  values  when  trying  to  approxi- 
mate the  desired  output  pulse. 
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A new  parameter  estimation  algorithm  may  need  to  be  devel- 
oped if  the  present  estimator  continues  to  fail  after  the 
above  recommendations  have  been  investigated.  Captain  J. 
Gary  Reid,  USAF,  has  a thesis  student  presently  working  on 
a function  minimizing  algorithm  which  could  be  of  use  in  a 
new  estimation  program.  For  more  information  he  can  be 
contacted  at  AFIT/ENG,  Wright-Patterson  AFB,  Ohio  45433 
or  AV  785-3450. 

As  noted,  none  of  the  test  cases  approximated  the  de- 
sired rectangular  pulse  very  well.  The  network  may  be 
changed  in  one  or  more  of  the  following  ways  in  an  attempt 
to  get  a more  rectangular  output  voltage  pulse: 

1.  Increase  the  number  of  meshes  in  the  network. 

2.  To  stop  postpulse  oscillations  in  the  load,  put 
a switch  between  the  network  and  the  load  such 
that  when  the  current  in  the  load  goes  negative 
the  switch  opens. 

3.  Distribute  the  stored  network  energy  throughout 
the  network  by  charging  all  the  capacitors  rather 
than  just  the  source  capacitors;  this  is  a con- 
ventional line-type  pulser  and  the  advantage  of  a 
lightweight  more  reliable  system  is  lost. 

4.  Put  a transformer  between  the  network  and  the  load 
and  have  the  turns  ratio  such  that  the  character- 
istic impedance  of  the  network  matches  the  load 
impedance;  this  could  decrease  the  pulse's  fall 
time  or  reduce  the  postpulse  oscillations. 

5.  Replace  the  E-line  with  a different  pulse-forming 
network . 
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Appendix  A Digital  Simulation  Program 


Operating  Instructions  and  Flow  Chart  for  RKFOUR 
and  Simulation  Program  Source  Listing 


Operating;  Instructions  for  RKF OUR 

Programmer:  B.  D.  Weathers 

Purpose:  RKFOUR  is  an  integration  routine  designed  to  provide  the  equivalent 

of  an  analog  computer  integrator  for  the  solution  of  normal -form 
differential  equations. 

Features 

An  arbitrary  number  of  integrators  can  be  used.  The  number  used  is  determined 

by  the  size  of  the  storage  arrays  set  in  the  calling  program. 

All  variables  that  affect  the  operation  of  RKFOUR  appear  in  the  calling 

statement;  multiple  uses  of  the  program  are  possible,  but  care  must  be 
used  to  prevent  time  desnehronization. 

The  independent  variable  step  size  is  determined  automatically,  but  if 

desired  this  program  can  be  made  to  run  with  a fixed  step  size. 

The  adaptive  step  size  feature,  if  used,  will  insure  little  trouble  with 

numerical  instability. 


The  limits  on  the  adaptive  step  size  may  be  set  in  the  calling  program, 

or  default  values  may  be  used. 


The  error  limits  that  determine  when  the  step  size  should  be  changed  are 

program  parameters,  and  hence  may  be  externally  controlled. 

An  adjustable  delay  is  available  to  prevent  eratic  step  size  changing. 

Once  the  criteria  is  satisfied  step  size  increases  by  the  appropriate 
factor  at  each  time. 

Step  sizes  may  be  changed  externally,  but  may  not  be  changed  in  the  middle 

of  an  iteration,  that  is  when  KEEP  is  equal  to  0. 

A solution  point  counter  is  contained  in  the  routine  to  enable  termination 
if  appropriate  independent  variable  criteria  is  not  available. 
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Suggestions  for  Use 

Before  using  this  subroutine  for  the  first  time,  a careful  examination  of 
the  flow  chart  should  be  made.  Also  a thorough  examination  of  theprogram  para- 
meter should  be  carried  out  to  insure  that,  the  default  values  of  the  parameters 
will  be  satisfactory,  if  used. 

You  must  set  all  non  zero  initial  integrator  outputs.  All  non-unity 
integrator  gains  must  also  be  set  in  section  2. 

TTje  value  for  H(l)  - Step  Size  - should  be  set  as  close  to  the  predicted 
sampling  interval  as  possible.  This  will  insure  a more  efficient  operation 
for  the  first  pass. 

If  running  in  "Open  Shop"  and  the  20  second  time  limit  is  just  a "bit  to 
short,"  add  "RUN=N0  CHECK"  to  the  control  card.  Do  this  only  after  your  name, 
after  you  are  SURE  the  program  is  free  of  programming  errors.  On  a large  job 
this  can  add  15%  to  the  available  computation  time,  by  eliminating  several 
compiler  tests. 


CALL  RKFOUR  (NINT,  Y,  X,  CA,  PO,  XX,  II,  ZZ1,  ZZ2,  ZZ3,  TIME,  PNTS,  KEEP,  IPL) 


NINT  - Number  of  integrators.  Storage  for  the  integrators  must  be  defined 
in  the  calling  program.  There  is  no  maximum,  but  the  example 
program  is  set  for  up  to  20. 

Y - An  array  that  contains  the  integrator  output.  Default  initial  values 
are  zero. 

X - An  array  that  contains  the  integrator  inputs. 

GA  - An  array  that  contains  constant  gains  associated  with  each  integrator. 
Default  values  are  1.0. 

PE  - A working  array  used  to  store  absolute  peak  values  of  the  integrator 
outputs.  This  array  is  set  to  zero  during  initial  entry  to  RKFOUR. 

XX  - An  array  used  in  two  ways.  If  I denotes  the  Ith  integrator,  then 
XX (1,1)  stores  the  past  value  of  the  integrator  output  and  XX (2, 1) 
stores  the  past  value  of  the  integrator  input.  During  the  integration 
cycle  the  XX(3,I)  and  XX (4,1)  positions  are  used  to  store  the  deriva- 
tive values  needed  by  the  integration  algarithm.  At  the  end  of  a 
cycle,  XX(3,I)  holds  the  new  integrator  output  and  XX(4,I)  holds  the 
new  integrator  input.  It  is  possible,  with  the  values,  to  implement 
a third  order  interpretation  for  the  values  of  the  integrator  outputs 
between  steps.  This  would  serve  to  offset  the  problem  of  obtaining 
values  at  discrete  time  points  and  would  be  more  efficient  than  step 
size  manipulation. 

H Array  - The  integration  control  parameters. 

H (1)  - The  current  independent  variable  step  size. 

H (2)  - Minimum  permissible  step  size  (Default  value  is  0.0). 

H (3)  - Maximum  step  size  (Default  value  is  1.0E  + 8). 
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H (4)  - Factor  used  to  modify  the  step  sizc-Must  be  greater  than  1.0 


(Default  value  is  SORT  of  2=1.414214). 

*» 

H (5)  - Proportional  to  the  lower  error  bound  (Default  value  is 
1.0E  - 4). 

H (6)  - Proportional  upper  error  bound  (Default  value  is  1.0E  - 02) 

H (7)  - Adaptive  step  size  control  (Set  to  -1  for  constant  step  size 
operation)  (Default  value  of  0.0  cause  adaptive  step  size 
operation) . 

H (8)  - Delay  in  increasing  the  step  size  after  error  criteria  is  met. 
(Default  value  is  3.0). 

H (9)  - Last  value  of  the  independent  variable. 

ZZ1,  ZZ2,  ZZ3  - Working  arrays  used  to  store  the  constants  needed  for 
the  integration  routine.  These  arrays  are 
TIME  - The  independent  variable  set  within  RKFOUR.  The  default  initial 
value  is  zero. 

PNTS  - The  counter  of  permanent  solution  extensions.  Use  this 

variable  for  termination  when  other  termiation  options  fail. 

KEEP  - This  variable  is  set  to  ONE  only  when  the  data  in  the  deriva- 
tive evaluation  section  is  permanent.  It  is  zero  at  all  other 
times.  This  variable  should  be  used  as  a control  upon  memory 
dependent  functions. 

IFL  - The  IFLag  array  (See  flow  chart) 

IFL(l)  - Set  to  -1  at  the  first  entry  to  RKFOUR,  and  is  0 at  the  first 
exit.  On  subseqeunt  exits  it  has  the  value  +1. 

IFL(2)  - Has  the  value  -1  to  +3  depending  upon  where  the  integration 
routine  is  in  its  operation  cycle. 
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This  flat  is  used  to  control  the  step  size  from  the  deriva- 
tive evaluation  section  of  the  calling  program.  It  is  set 
to  -1  for  step  size  reduction,  and  +1  for  an  increase  in  the 
step  size. 

This  flag  is  -1  when  initialializing  calculations  are  to  be 
made.  It  is  0 when  only  the  derivative  evaluation  is 
required,  and  is  +1  when  the  solution  has  been  successfully 
extended  over  one  step  of  the  independent  variable. 

This  flag  is  set  to  +1  immediately  after  IFL(4)  is  -1  for 
the  first  time.  It  is  intended  as  a signal  for  other  routines 
that  may  be  used  that  initialization  is  underway. 

This  flag  is  used  for  the  temporary  storage  of  the  variable 
KEEP. 
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The  Basic  RKI'OUR  Program  . 


DIMENSION  Y(20)4,  X(20)4>  GA(20)4>  PE(20) ,XX(4,20) ,H(9) 
DIMENSION  ZZ1  (20)  , ZZ2  (20)  , ZZ3(20)  , IFL(7) 


If  Other  Dimensioned  Variables  are  required  the  necessary  DIMENSION 
statements  go  here. 

SECTION  1 INSERTIONS  GO  HERE 

NRUN=  (Number  of  runs  to  be  made) 

DO  6 IRUN=1 ,NRUN 
IST0P=0 
KEEP=1 
IFL(l) =-l 

NINT=  (Number  of  integrators  required) 

1 IFL(2)=-1 

2 CALL  RKFOUR(NINT,Y,X,GA,PE,XX,H, ZZ1 , ZZ2, ZZ3, TIME, PNTS, KEEP, IFL) 
IF  (IFL[4])3,4,5 

3 INT=1 
TIME1=-1. 256789 
NC=1 

NCC=1  (Portion  of  Points  released  to  Section  5) 

FINTIM=  (Maximum  Value  for  TINE  in  this  run) 

PTMAX=  (Maximum  Number  of  solution  points  for  this  run) 

V, 

SECTION  2 INSERTIONS  GO  HERE 

4 IF (TINE  - TIME1) 7,8,9 

7 TIME1=TIME 

SECTION  3 INSERTIONS  GO  HERE 

8 CONTINUE 

SECTION  4 INSERTIONS  GO  HERE  (Exclude  everyting  nonessential) 

GO  TO  2 

5 IPR=0 
NC=NC-1 

IF  (NC)9,9,10 

9 NC=NCC 
IPR=1 

10  IF  (TINE . GE . FINTINI.  OR . PNTS . GE . PTMAX)  IST0P=1 
IF  ([IPR+ISTOP] .EQ.O)  GO  TO  11 

SECTION  5 INSERTIONS  GO  HERE 

11  INIT=0 

IF  (ISTOP) 1,1,6 

6 CONTINUE 
WRITE  (6,100) 

100  FORMAT  (1H1,  IX) 

STOP 

END 


Basic  RKFOUR  System 


There  are  a number  of  variables  in  the  RKFOUR  System  whose  value  depends 
upon  the  problem  undertaken. 

NRUN  - The  number  of  runs  to  be  made. 

NINT  - The  number  of  integrators  to  be  used  in  the  run  (limited  to  20  by 
storage  allocation  in  the  stock  system). 

NCC  - The  portion  of  solution  points  that  are  released  to  the  output 

section  (Section  5)  NCC  = 1 delivers  every  point,  NCC  = 2,  every 
second  point,  etc. 

FINTIM  - The  maximum  value  that  the  independent  variable  will  attain. 

The  program  will  terminate  upon  the  first  solution  point  to  meet 
or  exceed  this  value. 

PTMAX  - The  maximum  number  of  solution  points  in  any  one  run. 

IPNT  - The  point  counter  used  in  plot  routines.  It  may  be  set  to  0 in 
one  of  two  locations.  If  summary  curves  are  desired,  set  to  0, 
before  the  IRUN  DO  statement.  If  no  summary  curves  are 
desired,  set  to  0 after  the  IRUN  DO  statement. 

INIT  - This  variable  is  1 the  first  interation  step,  and  0 thereafter. 

This  can  be  used  as  a memory  flag  for  memory  dependent  functions. 

The  Basic  RKFOUR  System  contains  five  sections  where  the  particular  problem 
statements  may  be  inserted. 

1.  This  section  is  used  to  initialize  those  values  which  arc  not  run 
dependent,  or  whose  calculation  in  block  form  is  more  efficient.  In 

. this  case  they  may  be  indexed  by  IRUN  for  use  in  the  rest  of  the 
program. 

2.  This  section  is  used  to: 

a.  Set  non-zero  integrators  initial  conditions 
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and  are  evaluated  in  section  1,2,  or  3 as  is  appropriate. 

5.  This  section  is  used  to  call  the  output  routines.  If  NCC  is  not  equal 
to  1 it  is  not  entered  on  every  pass,  and  this  should  be  considered 
when  structuring  the  problem. 
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RKFOUR  Example 

Consider  the  simple  circuit  shown  below 


From  elementary  circuit  theory  we  know: 

C*VC  = 1R1  and  VR1  = V = VC 

Now  let  us  define  the  capacitor  voltage  as  a state  variable,  which  leads 
directly  to: 

Y (1)  - VC 

Therefore,  with  GA(1)  = 1.0/C 

VC=Y (1) 

1R1=(V-VC)/R 

X(1)=1R1 

are  the  derivative  evaluation  equations. 

In  the  example  program  on  the  following  pages,  the  capacitor  was  originally 
charged  to  -1  volt  before  the  switch  was  closed. 
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DIMENSION  Y(20) , X(20),  GA(20),  rr.(20) , XX(4,20),  11(9) 

DIMENSION  ZZ1(20),ZZ2(20) , ZZ3(20) , IFL(7) 

DIMENSION  AR(SOO) 

LOGICAL* 1 LR(400) 

REAL  IR1 

C SECTION  1 INSERTIONS  ******************************** 

C=I .0E-06 
R=1 .0E+04 

C IPNT=0  HERE  IE  SUMMARY  CURVES  ARE  DESIRED 

C NRUN= 

NRUN=1 

DO  6 IRUN=1 ,NRUN 
IPNT=  0 
WRITE (6 ,100) 

ISTOP=0 
KEEP=1 
IFL(l)  = -1 

C N1NT=  NUMBER  OP  INTEGRATORS 

NINT=1 

1 IFL(2)  = -1 

2 CALL  RKFOUR(NINT,Y,X,GA,PE,XX,H, ZZ1 ,ZZ2,ZZ3, TIME, PNTS, KEEP, IFL) 
IF(IFL(4))3,4,^ 

3 INIT=1 
TIME1=-1. 256789 
NC=1 

C NCC=NUMBER  OF  POINTS  KEPT 

NCC=1 

C FINTIM=FINISH  TINE 

FINTIM=0 . 5 

C PTMAX=N1AXIMUM  NUMBER  OF  POINTS 

PTMAX=400 

C SECTION  2 INSERTIONS  ****************************** 

Y(l)=-1 .0 
GA(1)=1 .0/C 
WRITE (6, 101) 

4 IF (TINE-TIME1)  7,8,7 

7 TIME1=TIME 

C SECTION  3 INSERTIONS  ****************************** 

V=1.0 

8 CONTINUE 

C SECTION  4 INSERTIONS  ************************ 

VC=Y(1) 

IR1=(V-VC)/R 
X(1)=IR1 
GO  TO  2 

5 IPR=0 
NC=NC-1 

IF  (NC)9,9, 10 

9 NC=NCC 
IPR=1 

10  IF (TIME . GE . FINTIM. OR . PNTS , GE . PTMAX) ISTOP= 1 
IF((IPR+ISTOP) .EQ.O)  GO  TO  11 

C SECTION  5 INSERTIONS  ******************************* 
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\ 

WR I TE ( 6 , 1 0 2 ) T I MR , I R 1 , VC 

44 

CALL  FP LOT ( 800 , 1 PNT , AR , LR , 0 , 1 , 2 , TIMR , V) 

45 

CALL  FP LOT (800,1  PNT , AR , LR , 1 STOP ,2,2, TIMR , VC) 

46 

11 

IN1T=0 

47 

IF(ISTOP) 1,1,6 

48 

6 

CONTINUE 

49 

WRITE (6, 100) 

50 

100 

FORMAT  (1111,  IX) 

51 

101 

FORMAT ( 5X, *TI ME * ,T25, 'RESISTOR  CURRENT* ,T50, 'CAPACITOR  VOLTAGE) 

102 

FORMAT(2X, 1PE12.5,T26,E1S.7,T52,E15.7) 

53 

STOP 

54 

END 

/DATA 

55 

SUBROUTINE  RKFOUR(N, Y,X, GAIN, PEAK, XX, II, Cl ,C2,C3, TIME, PUTS, KEEP, FL) 

56 

DIMENSION  Y (N) , X(N) , GAIN (N) , PEAK(N) , XX(4,N)H(9) ,C1(N) ,C2(N) , 

$C3(N) ,A(8) 

57 

INTEGER  F L ( 7 ) 

58 

DATA  A(l)/0.1/,A(2)/0.0/,A(3)/1.0E8/,A(4)/l. 414214/, A(S)/1. 03-4/ 

59 

DATA  A(6)/1.0E-2/,A(7)/0.0/,A(8)/3.0/ 

60 

IF (FL(1) ) 1,3,20 

61 

1 

DO  2 1=1, N 

62 

Y(I)=0.0 

63 

PEAK(I)=0.0 

64 

2 

GAIN(I)=1 .0 

65 

TIME=0 .0 

66 

PNTS=0 . 0 

67 

FL(4)=-1 

68 

FL(5) =1 

69 

FL(1) =0 

70 

DO  9 1=1,8 

71 

9 

H(I)=A(I) 

72 

RETURN 

73 

3 

Q=H(1) 

74 

FL(4)=1 

75 

FL(7)=H(8) 

76 

4 

H(1)=Q 

77 

FL(3)=0 

78 

DO  5 1=1, N 

79 

C2(1)=GAIN(I)*Q 

80 

C1(I)=0.5*C2(I) 

81 

5 

C3(I)=C2(I)/6. 

82 

IF (FL(1 ) )6 ,6, 10 

83 

10 

IF (FL(2) ) 31 , 33, 33 

84 

6 

FL(1)=1 

85 

7 

PNTS=PNTS+1 . 

86 

DO  8 1=1, N 

87 

XX(3,I)=Y(I) 

88 

8 

XX (4 , I)=X(I) 

89 

FL(4)=1 

90 

H(9)=TIME 

91 

RETURN 

92 

20 

IF(FL(3))25, 30,25 

93 

25 

IF(FL(2) -3)26,60,60 

94 

26 

IF(FL(3))21 ,30,22 

95 

21 

Q=H(1)/II(4) 
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IF(Q-II(2))23,4,4 

23  Q=ll(2) 

. IF(Q-ll(2))23f4,4 
22  Q=1I(1)  *11(4) 

IF(Q-I1(3))4,4,24 

24  Q=ll(3) 
IF(Q-II(]))4,30,4 

30  FL(5) =0 
IF(FL(2)) 31 ,40,45 

31  FL(6)=KFI3P 


• ^ • •**  * *,  ft. » 


o 


o 

0. 

2.C000300E-04 

-1.  D0GJ03GO00 

1.03000E-03 

1. 83 9 6750 E- 04 

-8.0967500E-01 

2.0G000E-03 

1.637  4 618E-Q4 

-6. 37^61806-01 

o 

3. 000  30E-0 3 

1.4316358E-04 

-*♦.61636841-01 

4.41421K-03 

1.2362-+43E-04 

-2.862-+431E-01 
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1.1 1 66 194E-04 
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8 . 4152797  E-05 
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2.631 37E-02 
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8.56037  34E-0 1 

O 
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1.17  367  01E-05 

8 , 82 13299E-01 

3 . 0 31 37E-Q  2 

9.650 1 639E-06 

9.0  349836E-01 

3. 31422E-G  2 

7.2728704E-06 

9. 272713GE-Q1 

<ff> 

3 ,597  96E-0  2 

5 • 43 1 2 17  IE-06 

9.45 187  63E-01 

3.87990E-02 

4. 1309331E-06 

9,536  3G67E-3 1 

4.16274E-02 

3.1132831E-06 

9.63867 1 2E-0 1 

Cj 

4.44559E-02 

2.3*+53374E-06 

9.7653663E-01 

4.E4559E-02 

1.5729338E-Q6 

9. 8h  27  0 16E-01 

5.24559E-02 

1. 054527  3E-06 

9, 694547  2E-01 

c 

5.64559E-02 

7,3695  512E-07 

9. 929304  5E -01 

6.04559E-02 

4.73942*+3E-07 

9.9526C58E-01 

6.44559E-02 

3. 1773038E-07 

9.9682269E-01 

o 

7.01128E-02 

1.806G163E-07 

9.93 1 9398E-0 1 

7.81128E-32 

8.1533633E-Q8 

9.9913416E-01 

3.61123E-02 

3.6354017E-08 

9. 9963146E-31 

9.41128E-02 

1.6648164E-08 

9, 9983352E-01 

1.02113E-01 

7. 52051 9 8 E- 09 

9. 9992479E-01 

. 

1.13427E-01 

2.52  i4G55E-09 

9. 9997477E-01 

«o 

1.29427L-Q1 

6 • 82  328  3 7E-1 0 

9 . 99993 1 8E-0 1 

1.52054E-01 

3.1295426E-10 

9. 9999687E-01 

1.840546-01 

5.72C0431E-10 

9.9999428E-01 

<r> 

2.29309E-01 

5.C315725E-09 

9. 9994993E-01 

2.74564E-01 

4.3733422E-08 

9, 9956267E-01 

2.97192E-01 

2.00  58598E-08 

9.9979941E-01 

vV 

3.19819E-01 

9.1999974E-09 

9. 999080  0E-01 

3.42-»47E-01 

4.219o344E-09 

9. 99957606-01 

3.650  74C-01 

1.9353603E-G9 

9.9993G65E-G1 

<v 

3.87702E-01 

3.8766493E-10 

9.9999112F-01 

4.10329E-01 

4.0713233E-10 

9 , 999  959  3E-0 1 

4.32957F-01 

1 , 367  3 396E-10 

9.99993131.-01 

t* 

4. 53534E-01 

3. 5646o5flE-ll 

9.9999914E-01 

4. 875  34E-01 

1.5654141E-10 

9.9599843E-01 

“( 

5.32839E-01 

1.3637  37  8E-09 

9. 9993631E-01 

<> 
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Flow  Chart  for  RKFOUR 
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Source  Listing  of  the  Simulation  Program 

These  job  cards 
are  needed  for 
the  digital 
simulator 
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PROGRAM  LPM(INPUT, OUTPUT) 

C 

C 

0»*»*»**»**»»V  »»*♦*■***»****»♦*♦»************■**»*♦**¥*♦. ♦♦♦♦♦.♦.TV.**** 

£*»*** ».♦».»».*».*..*♦.»***♦»*».*.♦*»»»**»*»**»**»*♦♦**.*♦»♦♦♦ ******** 
Q * * * » <•*♦♦«>**♦***¥»»*¥»**•»***♦¥*********.»*♦***•»*******♦**■**»****♦*♦#*** 
£**»*»**.*»»♦¥»*¥*¥♦***■»****.♦*♦♦»**»***.•♦****♦****»****»»*«**♦¥¥**♦♦* 

c 

c 

c this  program  simulates  th:  oischarge  OF  THE  E-LINi  PULSE-FORMING 
C NETWORK.  THIS  IS  DONE  3Y  FIRST  FORMING  AN  ANALOG  MODEL  OF  THE 
C NETWORK.  THE  MODEL  IS  THEN  USED  IN  A COMPUTER  SUBROUTINE  THAT 
C APPROXIMATES  AN  ANALOG  COMPUTER. 

C AN  ANALOG  MODEL  OF  THE  E-LINE  IS  DEVELOPED  IN  CHAPTER  IV  OF  THIS 
C THESIS.  THE  COMPUTER  SUBROUTINE  RKFOUR  APPROXIMATES  AN  ANALOG 
C COMPUTER.  THEREFORE,  IF  THE  EQUASIONS  USED  TO  DEVELOP  THE  ANALOG 
C MODEL  ARE  PROGRAMMED  INTO  RKFOUR  THE  DISCHARGING  E-LINE  PFN  CAN  3E 
C APPROXIMATED  (ScE  THE  OPERATING  INSTRUCTIONS  FOR  RKFOUR). 

C 

C THE  MAIN  PROGRAM  WHICH  CALLS  RKFOUR  HAS  SEEN  SET  UP  TO  HANDLE 
C 1.  VARYING  PULSE  RATES 

C 2.  VARYING  NUMBERS  OF  MESHES 

C 3.  VARYING  INITIAL  VOLTAGES  ON  THE  SOURCE  CAPACITOR 
C THE  CALLING  PROGRAM  HAS  ALSO  OECN  PROGRAMED  TO  MODEL  THE  SWITCHING 
C THAT  TAKES  PLACE  IN  THE  NETWORK  AS  THE  NETWORK  DISCHARGES.  RECALL, 
C A SWITCH  3ETWEEN  CS  AND  L ( 1 > CLOSES  AT  THE  START  OF  THE  PULSE, 

C OPENS  WHEN  THE  CURRENT  IN  L(l)  GOES  NEGATIVE,  CLOSES  AGAIN  WHEN 
C VCS  HECOMES  LARGER  THAN  VCi,  AND  OPENS  AGAIN  WHEN  THE  CURRENT  IN 
C L ( 1 ) AGAIN  GOES  NEGATIVE.  THIS  SWITCHING  CONTINUES  UNTIL  VCS 
C REACHES  A PREDETERMINED  VALUE  (SEE  PAGES  26-23  OF  THIS  THESIS). 

C 

C THE  CALLING  PROGRAM  CONTAINS  FIVE  SECTIONS,  A PROBLEM  STATEMENT  IS 
C PLACED  IN  A PARTICULAR  SECTION  DEPENDING  ON  THE  PURPOSE  OF  THE 
C STATEMENT.  SEE  THE  OPERATING  INSTRUCTIONS  FOR  AN  EXPLANATION  OF 
C EACH  SECTION. 

C 

C 

c**..**..  ♦♦***»»**..♦♦■,»*♦****»********»****•»*»****#*****♦♦**»******** 

(J*.****»».*****'.****»******4*************»*************  *************** 

c 

c 


LWV,T25,  107  0 , CM60CD 0 , STCSU, T7304O3, V ANNOY, 63835 
FTN, PL =8 000, OPT =3, I~ INPU T , L=OU T PUT , 3=LG0. 

LGO. 


OIMENS I ON  Y ( 20  ) , X ( 20 ) , GA ( 20 ) ,PE(2C) , XX (4 , 20 ) , M ( 9) 
DIMENSION  C (20) ,R(20) ,MFL (5) 

OIMEIlSION  271(20)  ,ZZ2(  20)  ,ZZ 3(20)  ,IFL  (7) 

REAL  IL1,IL2,IL3,IC1,IC2,IC3,IRL,L(20' 


MANY  OF  THE  VARIABLES  IN  THE  PROGRAM  ARE  OEFINEO  IN  THE 
OPERATING  INSTRUCTIONS  FOR  RKFOUR.  THOSE  VARIABLES  WHICH  HAVE 
HOT  BEEN  OEFINEO  ARE  OF.FINEO  OELOWI 

RL-- T ML  RESISTANCE  OF  THE  LOAD  RtSISTOR  IN  OHMS. 
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C L ( J) THE  INDUCTANCE  OF  THE  INOUCTOR  IN  THE  J'TH  MESH. 

C (HENRYS) 

C J=  1, 2, 3 ,••  • M J WHERE,  M=THE  NUH9ER  OF  MESHES. 

C v 

C . . C(J) THE  CAPACITANCE  OF  THE  CAPACITOR  IN  THE  J'TH  MESH. 

C (FARADS) 

C 

C R(J) THE  RESISTANCE  OF  THE  RESISTOR  IN  THE  J'TH  MESH. 

C (OHMS) 

C 

C PULSEC— THE  NUMBER  OF  PULSES  PER  SECOND 

C 

C TSTEP THE  TIME  BETWEEN  THE  BEGINNING  OF  ONE  PULSE  ANO  THE 

C BEGINNING  OF  THE  NEXT.  (SECONDS) 

C =1/PULSSC 

C 

C NFL(K)—  USED  IN  THE  PREPULSE  TIME  STEPPER  (SEE  SECTION  3) 

C K = l»  2, 3 

C 

C TIMELIO-USFO  IN  THE  PREPULSE  TIME  STEPPER  (SEE  SECTION  3) 

C 

C TIMES USED  IN  THE  PREPULSE  TIME  STEPPER  (SEE  SECTION  3) 

C 

C XMIN USEO  IN  THE  PREPULSE  TIME  STEPPER  (SEE  SECTION  3) 

C 

C ITRIG USEO  IN  SECTION  4 WHEN  MODELING  THE  SWITCH  IN  THE 

C NETWORK 

C 

C ZVCS USEO  IN  SECTION  4 WHEN  MODELING  THE  SWITCH  IN  THE 

C NETWORK 

C 

g*********»**¥¥**¥****** 
Q*.¥**¥**¥**4HHMM»»**¥***¥**¥*****»**«MKM,*¥******¥  *•*♦•¥**¥*¥*******¥*•¥♦ 

c 

c 

c 

c 

g * » * * ♦ **♦♦¥*♦■**♦*♦***••♦ 

C *****  SECTION  1 INSERTIONS  ***** 

C THIS  SECTION  IS  ENTERED  ONCE  PER  PROGRAM. 

C*****#¥¥¥4**¥#¥¥***¥*** 

c 

c 

c parameter  estimator  l and  c values 

c 

RL  = .165 
CS  = 1.333E-C4 
L ( 1)  = 2. 150E-07 
L(2)  * 5.244E-98 
L(3>  = 1 • 516E-07 
L ( 4 ) = 1.639S-09 
L(5>  = 1.270E-07 
C(l)  = 7.S26F.-06 
C (2)  = 7. 779 L- 06 
C ( 3 > = 8.515E-06 
C(4)  = 1.179E-05 
C(5)  = 2.029E-05 
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115 


120 


125 


130 


135 


140 


145 


150 


R(l)  = 

R(2 ) = 

R ( 3 ) = 

R(4I  = 

R(5>  = 

C 

PULSEC  = 100.0 
TSTEP  = 1,0/PUlSEC 
NFL  (1>  = 0 
Nn  (2)  = 0 
NFL (3)  = 0 
T I MEL  1 9 = C.  0 
TIMES=TSTEP 
ITRIG  = 1 
ZVCS  = 40.CE-06 
C 

XMIN  = 5.0E-06 

C THE  NUMBER  OF  RUNS  - 

NRJN=i 

00  6 IRUN= 1 , NRUN 
IPNT=0 
PRINT  100 
1ST  OP=0 

keep=i 

IFL (1 ) =-l 

C NINT  IS  THE  NUMBER  OF  INTEGRATORS 
HINT  = 11 

1 IFL (2 ) =-l 
C 

C 

C CALLING  THE  SU9RO JTINE  RKFOUR 

2 CALL  RKF0UR(NINT,Y,X,GA,P2,XX,H,ZZ1, ZZ2 , ZZ3 , TIME , PNTS , KEEP , IFL > 
C 

IF(IFL(4) ) 3,4,5 

3 INIT=1 

TIMti  = -1.256789 
NC  = 1 


15.CE-03 
13.0E-03 
1:>.  UE-03 
15. 0C-03 
15.0E-03 


155 


C 

C 

C 

C 


NCC  IS  THE  NUMBER  OF  POINTS  KEPT 
HCC  = 1 

FIHTIM  IS  THE  FINISH  TIME 
FINTIM  = 0.5 

PTMAX  IS  THE  MAXIMUM  NUMBER  OF  POINTS 
PTMAX  = 2000 


c * * * * 

* * * 

160 

Q ***** 

SECTION 

C THIS 

SECTION  I 

c * * * * 

c 

* * * 

H(l)  = 

0.5E-06 

165 

CSIC  = 

30.PE+03 

Y ( 1 > = 

CSIC 

GA ( 1)  * 

-1. 0/CS 

GA12)  = 

1.0/U1) 

GA  { 3 ) = 

1.0/C(1> 

170 

GA  ( 4)  = 

1. 0/L ( 2) 

GA ( 5)  = 

1.0/C(2) 

2 INSERTIONS 
i ENTERED  ONCE  P 

***** 


***** 

R RUN. 

* * * 


******* 
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GA  (6)  = l.C/L(3) 

GA  ( 7)  = 1,0/CHI 
GA  ( 9)  = 1.1/LIM 
GA (9)  = 1 • 0/C (4) 

• GA ( 13 ) = 1. 0/1(5) 

GA ( 11)  = la  0/C (5) 
PRINT  131 

4 IF(TIHE-TIMEl)7f8»7 
7 TIM£t=TIME 


* * ♦ * * 


*********** 


*****  SECTION  3 INSERTIONS  ***** 

THIS  SECTION  IS  ENTERED  TWICE  PER  TIME  STEP. 

A PREPULSE  TIME  STEPPER  HAS  BEEN  PROGRAMED  HERE.  THE  PRZPULSE 
STEPPER  MAKES  SURE  THAT  THE  OERITIVES  ARE  ACCURATE  JUST  PRIOR 
TO  THE  NEXT  PULSE  AND  THAT  THE  OUTPUT  VOLTAGE  AND  CURRENT  ARE 
CALCULATED  EXACTLY  AT  THE  BEGINNING  OF  A NEW  PULSE.  THIS  IS 
DONE  OY  FORCING  THE  SUBROUTINE  TO  MADE  TWO  CALCULATIONS  JUST 
PRIOR  TO  A NEW  PULSE  AND  MAKE  ONE  CALCULLATION  AT  THE  TIHE  OF 
A NEW  PULSE. 

*********************** 

IF ( NFL ( 2 ) • £Q , 2)  GO  TO  16 
IF (NFL (1 ) • EQ • 2)  GO  TO  15 
IF(TIMES-TIME.LT.XHIN)  GO  TO  24 
GO  TO  25 

24  IF(TIM£.EQ. TIMES)  GO  TO  16 

TIME  = TIMES-ADS(TIMES-TIM'E) /2.0 
NFL ( 1)  = NFL  ( 1 ) + 1 
GO  TO  25 

15  TIME  = TIMZS-U.0E-08) 

NFL  (2)  = NFL  ( 2 ) ■*•  1 

GO  TO  25 

16  TIME  = TIMES 

NFL  (3)  = NFL (3)  + 1 

IF ( NFL ( 3 ) • LT • 2)  GO  TO  25 

NFL (1)  = 0 

NFL (2)  = 0 

NF L (3  ) = 0 

TIMES  - TIMES+TSTEP 

25  CONTINUE 
29  CONTINUE 

8 CONTINUE 


*****  SECTION  4 INSERTIONS  ***** 

THIS  SECTION  IS  tNTEREO  4 TIMES  PER  TIME  STEP. 

THs.  SWITCHING  THAT  OCCURS  IN  THE  NETWORK  IS  MOOELEO  IN  THIS 
SECTION. 

********************** 


VO-  • 


■Bounanm 


PROGRAM  LPN 


74/74  OPT  = 2 


F T N 4 • 6 ♦ 4 46 


IL4  = Y (8) 

VC<*  » Y(9> 

IL5  = Y ( 1 0 > 

VC 5 = Y C 11 » 

C . . 

C THIS  IS  THE  SWITCH  MOOEL 

c 

C 

TI1EL10.  = TIMES-TSTEP 
IF(TI1EL10*TIHE.LT.ZVCS)  GO  TO  45 
GO  TO  46 

45  ITRIG  = 1 
GO  TO  47 

46  ITRIG  = 0 

47  CONTINUE 

IFIIL1.LE.C.O.ANO.TIME.GT.TIMEL10+2.0E-06)  GO  TO  30 
GO  TO  31 

30  IF(VC1.GT. VCS.AND. ITRIG.  EQ.l)  GO  TO  32 
IF (ITRIG. ECU  0)  GO  TO  33 

31  VCS  = Y ( 1) 

IL1  = Y ( 2) 

X C 1 » = IL1 
GO  TO  34 

32  Y ( 2 ) = 0.0 
IL1  = Y (2) 

X C 1 » = IL1 
GO  TO  34 

33  Y ( 1 > = CSIC 
VCS  = 0.0 

Y < 2 » = 0.0 
IL1  =0.0 
X(l)  = ILi 

34  CONTINUE 
• C 

c 

VR1  = Y ( 2) #R (1) 

VL 1 = Y ( 1) - Y (3) -VR1 

X ( 2 ) = VL  1 

IC1  = Y (2)-Y  (4) 

X ( 3 ) = IC1 

VR2  = Y(4)*R<2) 

VL2  = Y ( 3) - Y (5 ) -VR2 

XC4 ) = VL2 

IC2  = Y ( 4) -Y (6)  . 

X (5)  = IC2 
VR3  = Y(b>*R(3) 

VL 3 = Y ( 5)  -Y  (7)  -V.R3 

X (6)  = VL  3 

ICJ  = Y <6 J -Y  (8) 

X ( 7 ) = IC3 
VR4  = Y ( 3 » *R (4) 

VL4  = Y (7) -Y (9) “VR4 

X ( 9 > = VL4 

IC4  = YO)-Y(IO) 

X(9)  = IC4 
VX5  = Y ( 10 ) *R (5  > 

VLO  = Y(9)-Y(ll)-V^5 


c 
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X(10)  = VL5 
ICS  = Y{10)-(Y(11>/RU 
X(li)  = IC5 
IR'  = Y(11)/RL 
GO  TO  2 
5 IPR=0 
NC  = NC-1 
IF(NC) 9,9,10 
9 NC  = NCC 
IPR=1  • 

10  IFUIMC.GE.FINTIM.OR.PNTS.GE.PTMAX) ISTOP=l 
IF( (IPR+ISTOP) ,EQ. 0)  GO  TO  11 


* * * * 


*»*•****♦ 


*****  SECTION  5 INSERTIONS  ***** 

THIS  SECTION  CONTAINS  THE  OUTPUT  QUANTITIES. 

****************** 


* * * 


PRINT  102, TIME, VCS,Il_l,VCi,VC5,IRL 
11  INIT=0 

IF (ISTOP) 1,1,6 
6 CONTINUE 

100  FORMAT  < 1H1) 

101  FORMAT (6 X,* TIME*, 1 3X , * VCS* , 14X , * IL 1 * , 14X , ♦ VC1* , 14X ,* V C5 * , 14X , 

$ + 1 RL* ) 

1 C 2 FORMAT  C2X,1PE12. 5, 5X,1P£ 12.5, 5X,1PE1 2.5, 5X,1PE 12.5, 5X,1PE12.5, 
S5X, 1PE12.5) 

STOP 

ENO 
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SUBROUTINE  RKFOUR (N,Y,X, 6A I N , PE AK , XX , H , C 1 , C2 , C 3 , T IME , P NTS , KtE P , FL > 
DIMENSION  Y(N) ,X (N> , GAIN ( N) ,PEAK(N) ,XX(4,N> , H ( 9 ) ,Cl(N) ,C2 (N) , 

2C3 ( M) , A < 3 » 

INTEGER  FL (7 ) 

OAT  A A(i)/0.1/,A(2>/0.0/,A(3)/i.0E3/, A (4) / 1 . 41 421 4/ , A ( 5 ) / 1. OE-4/ 

’ OATA  A(5)/1.0->2/,A(7)  /Q . 0/ , A ( 3) /3. 0/ 

IFlFL(t) >1,3,20 
00  2 1=1, N 
Y(I)=0. 0 
PEAK(I)=0.0 
GA I N ( I ) = 1 . 0 
TIME* 3*0 
PNTS=0. Q 
FL ( 4) =- 1 
FL (5) =1 
FL  ( 1)  = 0 
00  9 1-1,8 
H < I ) = A 1 1 » 

RETURN 
Q=H (1) 

FL ( 4)  = 1 
FL  (?)  = H ( 8) 

H(1)=Q 

FL ( 3) = 0 

DO  5 1=1, N 

C2 ( I) =GA IN ( I) 4Q 

C1(I>=0.5*C2(I> 

C3  ( I) =C2 ( I) /6. 

IF (FL(1) >6,6,10 
IF(FL(2) > 31,33,33 
FL ( 1 ) = 1 
PNTS=PNTS+1. 

00  8 1=1, H 
XX ( 3, I ) =Y ( I ) 

XX(4,I)=X(I) 

FH4)  = 1 
H(  9 > =T IME 
RETURN 

IF ( FL ( 3 ) >25,30,25 
Ir(FL(2> -3> 26,60,60 
IF  (FL (3))21,30,22 
Q=rt(l>/H(4> 

IF ( Q-H ( 2 > ) 23,4,4 
G=H(2) 

IF(Q-H(1>><*,30,4 
Q=H (1 > *H (4) 

IF ( O-H ( 3) ) 4, 4, 24 
Q=H  (3) 

IF(Q-U(1> > 4,30,4 
FL ( 5>  = 0 

IF (FL (2) >31,40,45 
FL (6) = KEEP 
K£2P=0 
DO  32  1=1, N 
XX ( 1 , 1 > =XX (3,1) 

XX (2, 1) =<X  (4,1) 

FL  ( 2)  = 0 
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C 


V 


i 


f 


k" 

r 


c? 


P 


60 


. 65 


70 


75 


80 


85 


90 


95 


100 


105 


TIHE=H(9) +0.5+ H(l) 

00  34  1 = 1,  N 

34  Y ( I > = XX  <1,I)+C1<I)*XX(2,I) 

35  KEEP=0 
36 . . FL<4) =0 

RETURN 

40  FL  ( 2)  = 1 

DO  41  I = 1 , N 
XX  ( 3, 1 ) -=X  ( I ) 

41  Y(I)=XX(1,I)+C1(I)+X(I) 

GO  TO  35 

45  IF (FL(2) -2) 46,43,60  * 

46  FL ( 2) = 2 
TI'1E=H  ( 9)  + H (1) 

OO  47  1 = 1,  N 
XX(4,I)=X(I) 

47  Y(I)=XX(1,I)+C2(I)*X(I) 

GO  TO  35 

43  FL ( 2)  = 3 

IF ( H( 7 ) ) 49 , 52, 52 

49  KE£P=FL ( 6) 

50  00  51  1 = 1, N 

51  Y(I)=XX(1,I)+C3(I)MXX(2»I)*X(I)+2.MXX(3»I)+XX(4,I))) 
GO  TO  36 

52  ERR=0 . 0 

DO  56  1=1, N 

EST  ERR  = C2  (I)+ABS(XX (2,1) + X ( I ) - 2 • + XX (4,1)  ) 

IF(PEAKU)  ) 53,  53,54 

54  EST ERR =E5T ERR/ PEAK (I) 

53  IF ( ERR- E ST ERR) 55,56,56 

55  £RR=ESTERR 

56  CONTINUE 

IF ( ERR-H  (6) ) 57, 57 , 66 

57  IF(ERR-H(5) >53,59,59 

58  FL(7)=FL(7)-1 

IF ( FL ( 7 ) ) 65,65,49 
65  FL ( 3) = 1 

GO  TO  49 

59  FL ( 7 > =H ( 3) 

GO  TO  49 

66  FL (7) =H ( 3) 

GO  TO  21 

60  IFCH17) ) 7,61,61 

61  DO  63  1 = 1,  N 
ERR=43S ( Y ( I ) ) 

IF ( ERR- PEAK  ( I > ) 63 , 63 , 62 

62  PE1K(I)=ERR 

63  CONTINUE 
GO  TO  7 
END 


SYMBOLIC  REFERENCE  HAP  (R=l) 


c 


c 


» 
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Typical  Output  from  the  Simulation  Program 


rue 

vcs 

0. 

3.09COGE+04 

2.209702-08 

2.  999972+04 

2.600332-08 

2. 9 93962  *•  04 

2.697 982-0  8 

2 ,999962404 

2.79564ir-08 

2. 9999624C4 

2.  893292-03 

2.993962404 

2. 990  962- 3 3 

2.999952404 

3.G  58602-08 

2.999962404 

3.18622  2-C8 

2.999952+C4 

3.283922-08 

2. 399942404 

3.381 572-98 

2.999942404 

3 • 479  232-0  8 

2.999942404 

3.576882-08 

2.993932404 

3.67454E-G8 

2. 999332404 

3.772192-08 

2.999932404 

3 . 8o9?52-9  8 

2.999922404 

3.967512-08 

2. 999922404 

4. 065162-06 

2.999912404 

4. 162322-08 

2. 999912*04 

4. 269472-08 

2.999912404 

4.358132-08 

2.999932404 

4.456232-08 

2.999892404 

4.634342-08 

2.999892404 

4, 772  452-0  8 

2.999832404 

4.910502-08 

2.999872404 

5.04866008 

2. 399872404 

5.186762-96 

2.939352494 

5 . 324  3 700  8 

2.999852404 

5.46298008 

2.999842404 

5.  60  1 3800  3 

2.  3996 -*24 04 

5.739192-08 

2.999832404 

5. 934  502- J 8 

2.999822404 

6.129 81 C-08 

2. 999602+94 

6.325122-08 

2. 399792+04 

6. 520442-08 

2.999782+04 

b. 715752-08 

2.99976-+04 

6.911 06r>0  8 

2. 999752+04 

7.106372-08 

2. 999742+04 

7.301682-08 

2.399722+04 

7.496992-0 8 

2.999712+04 

7.692302-08 

2.999692+04 

7.837622-08 

2.993632+04 

8.052332-08 

2 . '399082  + 04 

8.359142-03 

2.999642+04 

8.635352-08 

2. 999612+  04 

8. 911562-38 

2. 999592+04 

9.187782-08 

2.999562+04 

9.4o3392-08 

2. 399532+04 

9.740202-08 

2. 399502+04 

l.C01»42-07 

2. 399462+04 

1.G29242-C7 

2. 999452+04 

1 . 0563*2- J 7 

2. 999422  + 0-4 

1.044512-37 

2. 999392+04 

1* 112132-0  7 

2.919352+04 

1. 15119* -07 

2. 993312+04 

1 . l*»9 75t-07 

2.999262*04 

».7.’9U'-J7 

2. 999212+04 

1.7'  <<•  -i|7 

2.9931 5fc.fr 04 

1.  Hf**'*!/ 

2>  '99  )11 2+  0-» 

1 -*r 

2.  M9Gj2«04 

|. 4' -If 

2.  ) Wf  BEfrOb 

1.4  .*-87 

2.  191«l4t*f4 

. . . . *.  „ . * 

» ««.a«r*a* 

IL1 

VC1 

0. 

0. 

3.060772+03 

4. 350142+  00 

3. 624822*03 

6.023252+00 

3.763302*03 

6.433922+00 

3. 896772+03 

6.961542+00 

4.0  52732*03 

7 • 456122  + 00 

4.168672*03 

7.967642+00 

4.304612*03 

.8.496102+00 

4.440332*03 

9. 041512+00 

4. 57645F+03 

9.603862+00 

4.712352*03 

1.013312+01 

4.843232+03 

1.077942+01 

4.984112+03 

1.139252+01 

5.119972+03 

1.202262+01 

5.255322+03 

1.266962+01 

5.391662*03 

1.333352+01 

5.327492+03 

1.401442+01 

5. 6u3332*03 

1.471212+01 

5.799102*03 

1.542682+01 

5.934892*03 

1.615342+01 

6.070662*03 

1.690602+01 

6.262662+03 

1.799422+01 

6.454532*03 

1.911542+01 

6.646572+03 

2. 027032+01 

6.83348&+C3 

2.145912+01 

7.03037**03 

2.263162+01 

7.222232*03 

2.393732+01 

7.414062*03 

2. 52273c+01 

7.6C5872+03 

2.655152+01 

7.797652+03 

2.790902+01 

7 . 9894  02*03 

2. 930022+01 

8.260532*03 

3. 132512+01 

8.5316JE+03 

3. 341742+01 

3.802622+0 3 

3.557702+01 

9.073572+03 

3.780402+01 

9.344472+03 

4. 009822+01 

9.613312*03 

4.245962+01 

9.386092+03 

4.488822+01 

1.015682+04 

4,738402+01 

1 • 042752*04 

4.994692+01 

1.069812+04 

5.257692+01 

1.  09t>86E  + C-fr 

5.527392+01 

1.123912+04 

5. 803792+01 

1.162132+04 

6.206112+01 

1.203372+0-fr 

6.621B12+01 

1.23c  ; JE  + Ofr 

7.050832+01 

1.275792+34 

7 • 493302+01 

1.314972+04 

7.949072+01 

1. 35314E+C4 

8, 418162  + Cl 

1.391302+04 

8.9005/2+01 

1 . 4 2945F  + 04 

9. 396232+01 

1.4r,753E  + 04 

9.905232+01 

1,303702+04 

1.042762+02 

1. 5433 J2+04 

1.096312+02 

1.597662*04 

1.174332+02 

1 . 631492+04 

1.254952+02 

1.7G529t+04 

1.338232+02 

1.759052+04 

1.424142+02 

1.812792+04 

1.512692+02 

1.866492+04 

1.603862+02 

1.923172+04 

1.697662+02 

1.97J31E  + 0-. 

1.794072+02 

9 n 77i,  i r 

0 

0 

« . «a7B^c*n? 

0 

• 

• 

VC5  IKL 

0.  0. 

- 0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

Q.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  G. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

C.  0. 

0.  0. 

0.  0. 

0 • 0 • 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  G. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  c. 

0.  0. 

I).  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  0. 

0.  -0. 

0.  0. 

0.  0. 

0.  c. 

0.  c. 

0.  0. 

0.  0. 

n . n . 
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r 

< 


C' 

5. 915976-07 

2.975965*04 

8 . 9732+5+  34 

3. 73  5396+C3 

0. 

0. 

6.993196-07 

2.  375436+04 

9.067756+04 

3. 8077)6+03 

0. 

0. 

7.C71326-07 

2. 974906*04 

9.155925+04 

3.8826C6+03 

0. 

0. 

■c. 

7.149+46-07 

2. 374366+04 

9.245756  + C4 

3.  957326  + 03 

0. 

0, 

7 . 227576-0  7 

2.973815+04 

9. 334246+34 

4.  03  5446  + 03 

0. 

0. 

7. 305696-07 

2.  )7  1265  + 04 

9. 422396+04 

4.  1U94.56  + 03 

0. 

0. 

T'V 

7.303326-07 

2. 972714+04 

9.510196+04 

4.185816+03 

. 0. 

0. 

7.461946-37 

2.972156+04 

9.597656+04 

4. 262546+03 

3 . 20  8676- C 4 

1, 9+4656-0  3 

7.5+0076-0  7 

2.-971536  + 04 

9.68474 f+C4 

4.339626+03 

7. 059086-04 

4. 278236-03 

7.613 196-07 

2. 971016+04 

9.771496+04 

4.417046+03 

1. 090956-03 

6. 611316-03 

PI  P 

7.696326-07 

2.97044-1  + 04 

9.857935+ J4 

4.49+796+03 

1.475996-03 

8. 945386-03 

7.77+4+6-07 

2.969366+04 

9.94396E+C4 

4. 572866+03 

1.661036-03 

1. 127906-32 

rf-, 

7 . 852576-07 

2.969276+04 

1.002975+05 

4.651236+03 

2.246076-03 

1. 361256-02 

7.930 596-07 

2.  96668c  + 04 

1.011506+05 

4. 729936+03 

2.631116-03 

1. 594616-02 

8.008326-07 

2.963096+04 

1.020006+05 

4. 8088 ac + 03 

3.016156-03 

1.827976-02 

■f' . 

8.086946-07 

2.967495+04 

1.028466+05 

4.883096+03 

3. 401196-03 

2.061336-02 

8.1)0  756-  j 7 

2.96733 5 +C+ 

1.029966+05 

' 4.902126+93 

3 . 469266- 0 3 

2.102586-02 

8.107666-07 

2.967336+34 

1.031706+05 

4.903146+03 

3 . 50896c- 0 3 

2. 126346-02 

3.114566-07 

2,367275+04 

1.031456+35 

4.910166+03 

3. 577036-03 

2. 167506-02 

8. 121 47r-3  7 

2, 967226+04 

l.C  321  36+J5 

4.923196+03 

3 . 645 J96-0  3 

2.209156-32 

3.128376-07 

2. 367176+04 

1.032945+03 

4.930216+03 

3.7 13166- C 3 

2.250406-02 

p 

8. 138 142-37 

2.967096+04 

1.033936+35 

4.9491p£+03 

3. 809426-03 

2. 308746-02 

8.191  356-07 

2.  3 669cc  + 04 

1 . 0 35486+  05 

4.954226+03 

3.9+5556-03 

2. 3912+6-02 

8.171 +3r-07 

2.966036+04 

1.  )3758F  + u5 

4.974126+03 

4. 138076-03 

2.507926-02 

r" 

8. 199106-07 

2 . 96662E  + C 4 

1.040356+05 

5.002296+03 

4. 410346-03 

2.672936-02 

8. 238166-07 

2.966316+34 

1.044756+05 

5.042196+03 

4.795396-03 

2.906296-02 

8. 277236-07 

2.966035+04 

1 ■ 043935+05 

5.06 ’156+  0 3 

5.180426-03 

3.139656-02 

c 

8.316296-07 

2. 96p7  06  + 34 

1.053116+05 

5.122156+03 

5 . 56546  -03 

3. 373006-02 

S • 355  <56-07 

2. 365395+04 

1.057206+35 

5.  162246  + 03 

5.9605C6-C3 

3.606336-02 

8. 39++if-07 

2.  965C85+  04 

1.0ol43E+C5 

5. 202376  + 03 

6.3355+6-03 

3.839726-02 

c 

d.  4 33+86-07 

2.  954766  + 04 

1 • 065536  + 05 

5. 242566  + 03 

6.720586-03 

4. 07  3086-0  2 

8.472p+£-07 

2 . 9644H6+  04 

1.069725+05 

5 .282606  + 0 3 

7 » 105u2E-03 

4.30b446-02 

9.51160E-07 

2. 964 1 +6  + G4 

1,07  38  56  + 05 

5.323096+03 

7. 49066-:- 03 

4.539796-02 

Ir 

•.539226-37 

2.963916+04 

1. 076766+05 

5. 35 1626+03 

7 , 7b2926- 0 3 

4.70+606-02 

9.563036-37 

2. 963606+04 

1. 078226 +C5 

5 • 365896+03 

7 . 955786- 0 3 

4.821t>S£-02 

8. 5 36 946-0 7 

2. 963695+  04 

1.07967E+05 

5. 3o017 6+03 

3.159986-03 

4. 94  5++6-0  2 

r~ 

8.580.356-07 

2. 963596+04 

1.081136+33 

5. 394456+03 

8. 364176-03 

5. C 6 92 06- 02 

K. 

8. 594  46..-0  7 

2. 363476+04 

1. 082586+05 

5 • 4 037  + 6+  G3 

8.565376-03 

5. 192956-02 

8.613)96-07 
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3. 062C5E-05 

1. 

1 0 36  02*  04 

4 , 367  36E+  04 

1. 193442*04 

8. 3/7752*03 

5.077422*04 

3.087052-05 

1. 

096472*04 

4, 298952*04 

1.092642*04 

8 . 42  4 562  *03 

5.105792+04 

3. 112055-05 

1. 

037432*0- 

4.22*112+04 

1. 3935/C+04 

8.460032*03 

5.127292*04 

3. 129732-05 

1. 

081922*04 

4.163612*04 

1.087782*04 

8.476702*03 

5.137392*04 

3.147412-35 

1. 

3 76422*04 

4.111952*04 

1.034052*04 

8.485572*03 

5. 1h2772  * 0 4 

3. 165082-05 

1. 

371012*04 

4.064902*0-* 

1. 079512*04 

8.4  86232*  03 

5. 143172*04 

3. 182762-05 

1. 

365672+04 

3.998 132* 34 

1. 074232*04 

8.478502*03 

5.13»492*04 

3.203  442-05 

1. 

360402*04 

3.942472*04 

1.168132*04 

8.462492+03 

5. 1 2 8782*0  4 

3. 211122-0  5 

1. 

05 *212*04 

3.333452*04 

1. 061302+04 

3.438542*03 

5. 114272*04 

3.235805-35 

1. 

053092*04 

3. 33677E+04 

1.05  >6  )2*04 

8.407212*03 

5. 095232*04 

3.25347F-05 

1. 

045032*  O4 

3. 73811c* 34 

1. 045292*04 

3 • 3c9262  + 03 

5. 072232*04 

3.271152-05 

1. 

040042*04 

3.7431 12*0-* 

1.030132*04 

5.325602*03 

5, 04»822*0  4 
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1. 

333102*04 

3. 7C24  6L*  04 

1.0261/2*04 

8. 277272+03 

5.  C16632*  1-4 

3.31333005 

1. 

023212*04 

3.663372*04 

1. 0110.62  + 04 

8. 203032*03 

4. 971572*04 

3.349HO05 

1. 

013592*0-. 

3.602572*04 

9. 86o3  92  + 03 

8 . C >2142  + 0 3 

4.904332+04 

3.384 542-05 

1. 

0090b2*C4 

3.573652*04 

9. 66946E  *03 

7.981722*03 

4. 837412+04 

3.4C9542-05 

1. 

0123  92*0-. 

3.564682+04 

9. 534062*03 

7.907812*03 

4. 792612*04 

3. 434542-05 

3. 

957042*03 

3. 56227c ♦ 34 

9.422022*03 

7. 839442*03 

4,  751 162  *04 

3. 459  >400  5 

9. 

39022c*  03 

3. 563632+04 

9.337752*03 

7.777712*03 

4.713772*04 

3.464540  J5 

3. 

323372*03 

3.585312*04 

9. 282222  *03 

7.723052*03 

4.681642*04 

3.5C9542-03 

9. 

7565 jt*  03 

3. 5b4  09E  *0-* 

9. 253032*03 

7 . 67519c*  0 3 

4. 651632+04 

3.534p4E-05 

3. 

68971E*03 

3. 567  27E+94 

9. 245222*03 

7.633292*03 

4.62c24C*04 

J. 552222-05 

9. 

5 4,6>c*  0 3 

3. 547942*  04 

9.24)182*03 

7 . 6C  '>362*0 3 

4.610032*04 

3.569 39r-05 

9. 

3956JI+0J 

3,534322+04 

9.2583  !E  + 03 

7 .581562+03 

4,594832*04 

3.587>7r-05 

9. 

5436.82*03 

3.51617E+C-* 

9.270432*03 

7.567652*03 

4.580392*04 
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9.5C343E-05 
9.528432.-00 
9.55343E-05 
9.57843E-Q5 
9.6P343E-95 
9.62343E-05 
9.646112-05 
9.663792-05 
9.681462-J5 
9.6+9142-05 
9.71632L-Q5 
9.734502-05 
9.75217E-05 
9.769352-05 
9.79*652-05 
9.819952-05 
9.844851-05 
9. 969052-35 
9, 8943KE-05 
9.919852-35 
9. 9+485E-35 
9,969352-05 
9. 994  852-0  5 
1. 001992-04 
1, 004492-04 
1.006992-04 
1. 009492-04 
1.011992-04 
1.  C1449E-04 
1.01699t>04 
1 ,019492-04 
1.  021992-04 
1.024492-34 
-.026992-04 
1. 029492-04 
1.031992-04 
1.034492-04 
1.036252-04 
1 . 038  0 22-0  4 
1.039792-04 
1.041562-04 
1. 043  J2E-04 
1.045  J«»c-04 
1.046362-04 
1.048632-04 
1.050402-34 
1.052162-04 
1.054662-04 
1.053202-04 
1.361 73E-34 
1.085272-34 
1. 06734F-04 
1. 068O12-9  4 
1. 070572-04 
1. 072  342-34 
1.074112-04 
1.07S33c-04 
1.077642-34 
1.07941" -04 
1.081182-04 
1.032952-04 
1.063+5F-04 
1.C83 982-04 
1.092522-04 
1.095022-04 

1. 


2.0  8 49*2+03 
2.071522+03 
2.  05  8 1 92  + 03 
2.344962+03 
2.031012+03 
2.013762+03 
2. 309592+03 
2.  0 004ot  + 03 
1. 991372+03 
1. 992322+C 3 
1.973322+03 
1.964352+03 
1.955432+03 
1.946542+03 
1.934042+03 
1.921612+03 
1.909262+03 
1. 396972+03 
1.  88477::  + 03 
1.372632+03 
1. 8SC57E+03 
1.348562+03 
1.63667c+03 
1.324832+03 
1.313072+03 
1. 301402+03 
1.76 '9802+03 
1. 77825E+03 
1.766842+03 
1.755482+03 
1.7441 9c. + 03 
1.732992+03 
1.721052+03 
1.710802+03 
1.699812+03 
1.683892+03 
1.673042+03 
1. 670412+03 
1.662822+03 
1. 655252+03 
1.647722+03 
1.64922E+G3 
1.632752+03 
1. 625312  + 03 
1.617912+03 
1.610532+03 
1.  60  3192  + 03 
1.592372+03 
1.579332+03 
1.564022+03 
1.549802+03 
1.5427*2+03 
1.535725+03 
1.52872b +03 
1.521772+03 
1.514342+03 
1.50795E+03 
1.501092+03 
1.494262+03 
1.48747E+03 
1.430702+03 
1.47119F.  + 0J 
1,457842+03 
1.44462E+03 
1.433332+03 
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7.183152+03 
7.132872+03 
7 • 0825 JE+03 
7.33277E+03 
6.933532+03 
6.935342+03 
6.901922+03 
6. 6691 4E  + 0 3 
6. 837032+03 
6. 905612  + 03 
6.774872+03 
6.74479E+03 
6.715322+03 
6.666402+03 
6.046312+93 
6.605952+03 
6 , 5630  3E+03 
6. 529332+33 
6 , 4905 12+03 
6.451432+03 
6.411932+03 
6. 371812+ 33 
6. 331102+33 
6. 2o 9 612  + 03 
6.2479 9 E +03 
6.2C530E+03 
6.163392+03 
6. 120932+03 
6.07  876E  + 0 3 
6.0  3696E  + 0 3 
5.995752+93 
5.95533E+J3 
5.915712+03 
5.9773 7 E+33 
5, S3937E+03 
5.6C233E+03 
5.7665U  + 33 
5.74160E+03 
5. 716332+03 
5.692402+03 
5.8630  92  + 03 
5.64388E+C3 
5.619732+03 
5.595p62+03 
5.57135E+03 
5.54703E+03 
5.52237l+03 
5.487722+03 
5.437862+33 
5.387372+03 
5.33&47L+03 
5.313962+03 
5.265432+33 
5. 263072+03 
5.234732+03 
5.209672+03 
5.184762+03 
5. 169 12E+03 
5.135762+03 
5.111722+03 
5. 008012  + 03 
5,0r>506F.+03 
5. 00961c+03 
4.965332+03 
4.93471E+03 
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2.02037E+03 
2.007802+03 
1. 995052  + 03 
1.982112+03 
1 • 96399E  + 0 3 
1.9557*2+03 
1.  9463  If.  + P 3 
1.936892+03 
1.927452+03 
1.91804E+C3 
1. 903672+03 
1.09936E+03 
1.89019E+03 
1.881102+03 
1. 863472+03 
1.856132+03 
1. B44C92+03 
1.632332+03 
1.820662+03 
1.8096iE+03 
1.798622+03 
1.787742+03 
1.776972+03 
1.766242+U3 
1.755502+03 
1. 744712+03 
1 • 73385E  + 0 3 
1.722632+03 
1.711812+03 
1.700642+03 
1.68 J39E+03 
1.676082+03 
1 • 6667  6E  + 03 
1. 655472+  03 
1.644242+03 
1.033132+03 
1 • 62216E+  03 
1.614522+03 
1.606972+03 
1.599532+03 
1.5921 3E+03 
1. 584932+03 
1,577822+03 
1.570792+03 
1, 563852+03 
1 . 55698E+  C3 
1. 550192+03 
1.540672+03 
1. 3273h2+03 
1.514072+03 
1. 500772+03 
1. 49419F.+  03 
1.487392+03 
1.480662+03 
1.1.73902  + 03 
1.467122+03 
1.460312+03 
1.45349E+U3 
1.  ■.45662  + 03 
1. 439832+03 
1.4330 12 +03 
1.423412+C3 
1. 409952+03 
1. 396712+03 
1.387592+03 
, iT«i.r,  r.m 


1.568212+03 
1.567912+03 
1 . 5 77t>52  + 0 3 
1 . 5l,/16E  + 0 3 
1 • 0667bL  + 0 3 
1 , 54  6 362  + 03 
1,539032+03 
1.531742+03 
1.62*492+03 
1 . 51729c  + 03 
1.510162+03 
1. 503092+03 
1.496102+03 
1.489182+03 
1 • 4V’9'552  + C 3 
1 . 47  0072  + C 3 
1.460762+03 
1.451 582  + 0 3 
1.442532+03 
1.463682+C3 
1. 424712+03 
1 . 41589E  + 0 3 
1. 457092+03 
1.358292+03 
1 . 36948C  + 0 3 
1.38  0652  + 0 3 
1.371792+03 
1 . 3b290c  + 0 3 
1 . 354912+03 
1. 345102+03 
1. 32621E+C3 
1. 327342+03 
1.318522+03 
1.3C978E+C3 
1.301132+03 
1 . 2925o2  + 0 3 
1 . 28*142  + 03 
1.276252+03 
1.272432+03 
1 . 2C666E  + 0 3 
1.260952+03 
1.255292+03 
1.249692+33 
1.244132+03 
1 . 238622  + 0 3 
1.233142+03 
1.227692+03 
1.220322+03 
1.209212+03 
1 . 198*42  + 03 
1 , 187  652  + 03 
1 . It  2252  + 0 3 
1.176352+03 
1. 171452  + C 3 
1.166352+03 
1. 160662+03 
1 . 15527E  + C 3 
1.149,892  + 03 
1. 144532+03 
1.139192+03 
1 , 133381  + 0 3 
1 . 126422  + C 3 
1.115992+03 
1.  105722  + 03 
1.098582+03 
< n o . 5 •».:  * n 1 


9.636112+03 
9. 62  3t7£  + 0 3 
9, 56092Z  + 0 3 
9. 497932+03 
9.434892  + 0 3 
9.371692+03 
9. 327462+03 
9. 283282+03 
9.239352+03 
9.19571E+03 
9. 152472+03 
9.109652+03 
9.067252+03 
9. C25352+03 
8.966952+03 
6.909532+03 
8. 853082  + 0 3 
6.797482+03 
8.742642+03 
8.636392+03 
9.634612+03 
8.531132+03 
8.527602+03 
8.4  74402  + 0 3 
8.42110F.  + 0 3 
8. 367562+03 
8.313672+03 
8.262032+03 
8.206112+03 
8. 152132+03 
8. C9322E+03 
6.0*4432+03 
7.991C62+03 
7.939092+03 
7.  9dii64E  + Q3 
7.833822+03 
7.732632+03 
7 ■ 74t>98£+  0 3 
7. 711682+03 
7.676722+03 
7.642112+03 
7.607842+03 
7.573892+03 
7.5*0212+03 
7.506732+03 
7.473572+03 
7. 440532+03 
7.  39  + 04E+  0 3 
7. 328562+03 
7 • 263252 ♦ C 3 
7.197682+03 
7.169172+03 
7.132452+03 
7.099722+03 
7. C6658E+0 3 
7.034202+03 
7. 001632+03 
b. 95905E  + 0 3 
6. 9 3657E  + 0 3 
6. 90*202+0 3 
6.872002+03 
6. 626762  + 0 3 
6. 763562+03 
6.701362+03 
6. b5  3032  + 0 3 
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2.729462-04 
2.731225-34 
2.73J72F-04 
2.73726F-04 
2.74C7U-04 
2.743292-04 
2.745792-94 
2.746?9i:-P4 
2.750792-34 
2.753292-04 
2.755792-04 
2.767562-04 
2.759332-04 
2.761102-04 
2*  762862-9  4 
2.7C463E-04 
2.7664PE-34 
2.768172-04 
2.769-342-34 
2.771702-04 
2.773472-04 
2.775372-04 
2.77951L-04 
2.78304E-04 
2.75658E-04 
2.738352-04 
2.799112-04 
2.791892-04 
2.793652-04 
2.795-.22-04 
2.797 192-34 
’•7 98 95 2- 04 
• 833  722-54 
2. 832492-04 
2.804262-04 
2.8C676E-34 
2.810292-34 
2.313832-34 
2.816332-04 
2.818832-04 
2.821332-04 
2. 82333c -34 
2.826335-0  4 
2. 82  9832-04 
2. 830  39“ -G4 
2. 8 32  365-0  4 
2. 834132-04 
2.835302-04 
2. 837675-04 
2. 839  435-3  4 
2.841202-34 
2. 842972-04 
2.864742-04 
2. 646505-04 
2. 84909F-0  4 
2.882342-04 
2.856C81 -04 
2.85961004 
2.861 352-04 
2. 863 15E-04 
2.  864'I1l-04 
2.8  .6682-04 
2.868-.5E-04 
2.870  22C-04 
. 87 199F -3  4 
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2. 125682*01 
2.  116032+01 
2.102572*01 
2. 383622*01 
2.  364352*  01 
2. 051672*01 
2.338532*01 
2.325572*01 
2.012642*01 
1.999792*01 
1. 98702E*01 
1.978042*31 
1. 969102*01 
1.960202*01 
1.951342*01 
1.94252E*01 
1. 933742*01 
1.324992*01 
1.916292*01 
1. 90762E*01 
1.398992*01 
1. 836852*01 
1.869822*01 
1.952932*01 
1.836202*01 
1.327892*01 
1 • 81 962E*  01 
1.311392*01 
1.903192*01 
1.795042*01 
1.736922*01 
1.778842*01 
1.770802*01 
1.762792*01 
1. 754832*01 

1. 7  4 3b2E  + C 1 
1.727902*01 
1.712332*01 
1,701402*01 
1,6905^2*01 
1.679742*01 
1.569022*01 
1.653362*01 
1, 64777C*C1 
1.640332*01 
1. 63291E*01 
1.625532*01 
1.61  1182*01 
1.510072*01 
1.603582*01 
1 • 593  33c.*  01 
1.599122*01 
1,5319 32*01 
1. 574782*01 
t . “647  22*31 
1.550602*01 
1,536612*01 
1.522742*01 
1.515862*01 
1.539012*01 
1.502192*01 
1.495402*01 
1 . 48365E*  01 
1.481922*01 
1.475232*01 

I.  i.  «.  a c * f «.  n , 


7.258412+01 
7.224692*01 
7.17751F+01 
7.111982*01 
7.047902*01 
7.003372*01 
6.959362*01 
6.915742*31 
6. 87240E+01 
6. 829232*91 
6.786322+01 
6. 756122+01 

6.7  26  08I.+  0 1 
6.696212*01 
6 « 666  3 3 E + 01 
6.637C6E+01 
6.6Q780E+01 
6.578  742  + 01 
6.549542*01 
6.521072*01 
6 , 4924  0E+  0 1 
6. 451942  + 01 
6.394762*01 
6.337442*01 
6.279912*01 
6.251072+01 
6.222222+01 
6.193342+01 
6.164472*01 
6.135642+01 
6.106892+01 
6.078262*01 
6.049772*01 
5.021472*01 
5.9934 jE*01 
5. 954132+01 
5.899592*01 
5.84622E+01 
5.809132*01 
5.772482*01 

5.7  361d2+  0 1 
5.700132+01 
5.664272*01 
5.628552*01 
5.b03402*Cl 
5.576352+01 
5.553402+01 
5.528552*01 
5.503332*01 
5.479232+01 
5.454742*01 
5.430352*01 
5.406052*01 
5.3S1532+C1 
5, 347o52*01 
5.299412*01 
5.251162*01 
5.202812*01 
5.17853E+01 
5. 154322*01 
5.130002*01 
5.105662*01 
5.031292*01 
5. 056932*01 
5 • 032602*01 
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2.  057932*01 
2. 048552+01 
2.035202*01 
2. 016352+01 
1.997632+01 
1. 964672  + 01 
1. 971842*01 
1.959232*01 
1. 346742+01 
1.934382+01 
1.922062+01 
1.913332*01 
1.904652+01 
’l. 395972+01 
1.867312+01 
1.878672*01 
1.370082*01 
1. 361 5oE+  01 
1.853102+01 
1. 844722*01 
1. 336432+  01 
1. 324842+01 
1.303702*01 
1.792792+01 
1. 777062+01 
1.769232+01 
1.761402+01 
1, 75361E  + 0 1 
1. 745822+01 
1.733022*01 
1.730222*01 
1.722412+01 
1.  71459E+01 
1.70u76E+Cl 
1.693922+01 
1.687842+01 
1.672222*01 
1.656762+01 
1.640962*01 
1.635232+01 
1.62*732*01 
1.614442+01 
1.604192+01 
1.594002*01 
1. 53630E  + 01 
1.579642+01 
1.572522*01 
1.565402*01 
1. 553312+01 
1.551252*01 
1.544232*01 
1 • 537  26E*  0 1 
1.530342*01 
1.523432+01 
1.513862*01 
1.500442*01 
1.467212*01 
1.474152+01 
1.467662*01 
1.461232*01 
1.454332*01 
1.443442+01 
1. 442052*01 
1.430632*01 
1.429312+01 
« «,  9 ?a . r a a. 


1.630692+01 
1.623722*01 
1.613362*01 
1.599922+01 
1.506402*01 
1.574722+01 
1.564252*01 
1.553332*01 
1.54.4122*01 
1.534272*01 
1.524302*01 
1 • 516742*  0 1 
1.510062*01 
1.503392*01 
1.496122*01 
1.489152*01 
1.432132*01 
1.475212+01 
1 . 4b8242  *01 
1.461272*01 
1.454302+01 
1 • 445672+0 1 
1.433482*01 
1.421282*01 
1. 4G621L+  0 1 
1.401822+01 
1.396012*01 
1. 369312*01 
1. 333522+01 
1.377432+01 
1.371332+01 
1. 3o5232+Cl 
1, 359132*01 
1.363032+01 
1. 346332*01 
1.338312+01 
1.326112*01 
1.313912*01 
1.305062*01 
1.296662*01 
1. 288452*01 
1.280442+01 
1.271812+01 
1.263192*01 
1.257242*01 
1.251862*01 
1.2h6342*01 
1.240832+01 
1.235312*01 
1.230382*01 
1.224352*01 
1. 219a32*  01 
1. 214402+01 
1.209172*01 
1.2C 1782+01 
1.191322+01 
1.  180872  + 01 
1. 169332+01 
1 . 104312+0 1 
1. 15938E+C1 
1. 154152+81 
1.  148922+01 
1. 143702+01 
1.  1 33472  + C 1 
1. 183242+01 

, « iftn.r.n, 


9. 882952+01 
9.840712*01 
9.780972*01 
9.695492*01 
9.606462*01 
9.543762*01 
9.430292*01 
9. 41d062+01 
9. 358322*01 
9.298582+01 
9.236352+01 
9.192352+01 
9.151872*01 
9.10  3622+01 
9.067382*01 
9. 026142*01 
8.982902*01 
8.940662*01 
8.893412+01 
8.83O172+01 
8. 813932+01 
3. 761652*0 1 
6.637732+01 
6.6138 IE +01 
8.  53*-o02  + 0 1 
8.495882+01 
8.460682+01 
8.42372E+01 
8.384992+01 
8.343032+01 
8.311072*01 
8.274112+01 
3.237142*01 
8.200182*01 
8.163222*01 
8. 110952*01 
8.037022*01 
7.963102*01 
7. 93 95 8E *01 
7.858562*01 
7.608772+01 
7.760232*01 
7.707962*01 
7.656692+01 
7, 61961E*  0 1 
7.587052*01 
7.553602*01 
7.520162*01 
7,4do72£*01 
7.455042*01 
7.423662*01 
7.391672+01 
7.359992*01 
7.328312*01 
7.283512+31 
7.220142+01 
7. 1567o2+0l 
7.089692+01 
7.066462+01 
7.026532+01 
6.994852+01 
6.963172+01 
6.931492*01 
6.894802+01 
6.668122+01 
r.  « mi.i.r  *n  « 


J. 620 326-04 

2.23520E+00 

7.424746+03 

2. 152136+C3 

1.719356+00 

1. 041856+01 

3. 622096-34 

2.22537E+30 

7.400256+03 

2. 145696+00 

1 . 7t  8386  + 0 0 

l.C 35696+ 01 

3.62365E-04 

2.215576+00 

7.373236+00 

2. 139426+00 

1. 703076+00 

1.032176+01 

3.625626-04 

2. 205316+03 

7.343656+00 

2.133036+00 

1.700176+00 

1. 030416+01 

3.627  396-04 

2. 196136  + 00 

7.311956+00 

2. 126036+30 

1.691*66+00 

1.025136+01 

3.629 166-04 

2.19  6+26+33 

7.273616+00 

2. 118576+00 

1 .682756+03 

1. C19856+01 

3.630926-34 

2. 176796+00 

7.244166+00 

2. 110496+00 

1.674036+00 

1. C14575+0  1 

3. 632095-04 

2. 167216+00 

7.209116+00 

2.101836+00 

' 1 • 665826 +00 

1. 009296+01 

3.634465-34 

2.157676+00 

7.173926+00 

2.092796+00 

1 . 656616+0 0 

1.  00+Clc+01 

3.636236-04 

2.148186+00 

7.133966+00 

2.083336+00 

1.647906+00 

9. 987256+00 

3.638  736-0  4 

2. 134846+03 

7.090476+03  . 

2. 069596+00 

1.635576+00 

9. 912576+00 

3.641235-04 

2. 121566+30 

7.043446+00 

2. 055636+00 

1.623256+00 

9. 837906+00 

3.  6*3735-04 

2.108426+00 

o .998036+00 

2. 041816+00 

1 • 6 13936+  0 C 

9. 7 5 3226  + 0 0 

3.6+6236-34 

2.095336+03  . 

6.9542J6+00 

2. 023056+00 

1 • 59361E  + 0 0 

9.638556+00 

3. 649766-04 

2.376976+00 

6.894366+00 

2. 0C365E+  00 

1 • 5 8119E  + 00 

9.582946+00 

3.663306-04 

2.053766+03 

6. 838976+00 

1.983936+00 

1.569376+00 

9.512646+00 

3. 656836-04 

2.040686+33 

6.707426+00 

. 1. 968476+00 

1. 560S66+CC 

9.45 97 36+00 

3. 699336-04 

2.027996+00 

6.734136+00 

1.954216+00 

1. 554706  + 0 C 

9.422406+00 

3. 661 1GE-04 

2 • 0 19f  46  + 00 

6.732316+00 

1.944346+00 

1.548896+00 

9.337206+00 

3.662876-04 

2.010136+00 

6.711036+00 

1.93+286+00 

1 . 543176+00 

9. 33*396+00 

3.664646-04 

2.001246+00 

6 . 69 1156+00 

1.924476+00 

1.534376+00 

9.299196+00 

3.636416-04 

1.992366+00 

6.6721 J6+C0 

1.91509E+0C 

1.531466+00  . 

9.231596+00 

3.668176-04 

1.983556+00 

6.653586+00 

1.906066+00 

1.525656+00 

9.2+0396+00 

Appendix  B Parameter  Estimation  Program 


Operating  Instructions  for  ZXMIN  and  MEXP, 

Flow  Chart,  and  Estimation  Program  Source  Listing 


on  nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnno 


Operating  Instructions  for  ZXMIN  (Ref  5) 


SUBROUTINE  ZXMIN 


(FUNCT,  N, NSIG, MAXFN,  IOPT,  X,  II,  G,  F,  W,  IER) 


FUNCTION 

USAGE 

PARAMETERS  FUNCT 


N 

NSIG 


MAXFN 

IOPT 


X 


II 


G 

F 

W 


- A QUASI -NEWTON  ALGORITHM  FOR  FINDING  THE 

MINIMUM  OF  A FUNCTION  OF  N VARIABLES. 

- CALL  ZXMIN  (FUNCT, N.  NSIG,  MAXFN,  IOPT, X, II, G, F, 

K.IER) 

- A USER  SUPPLIED  SUBROUTINE  WHICH  CALCULATES 

THE  FUNCTION  F FOR  GIVEN  PARAMETER  VALUES 
X (1)  ,X(2),„,,X(N)  . 

THE  CALLING  SEQUENCE  HAS  THE  FOLLOWING  FORM 
CALL  FUNCT  (N,X,F) 

WHERE  X IS  A VECTOR  OF  LENGTH  N. 

FUNCT  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT 
'IN  THE  CALLING  PROGRAM.  FUNCT  MUST  NOT 
ALTER  THE  VALUES  OF  X(I) , 1=1, , , , ,N  OR  N. 

- THE  NUMBER  OF  PARAMETERS  (I.E..  THE  LENGTH 

OF  X)  (INPUT) 

- CONVERGENCE  CRITERION.  (INPUT).  THE  NUMBER 

OF  DIGITS  OF  ACCURACY  REQUIRED  IN  THE 
PARAMETER  ESTIMATES. 

THIS  CONVERGENCE  CONDITION  IS  SATISIFIED  IF 
ON  TWO  SUCCESSIVE  ITERATIONS.  THE  PARAMETER 
ESTIMATES  (I,E,,X(I),  1=1,,,,, N)  AGREE. 

COMPONENT  BY  COMPONENT,  TO  NSIG  DIGITS. 

- MAXIMUM  NUMBER  OF  FUNCTION  EVALUATIONS  (I.E.. 

CALLS  TO  SUBROUTINE  FUNCT)  ALLOWED,  (INPUT) 

- INPUT  OPTIONS  SELECTOR . 

IOPT  = 0 CAUSES  ZXMIN  TO  INITIALIZE  THE 
HESSIAN  MATRIX  II  TO  THE  IDENTITY  MATRIX. 

IOPT  = 1 INDICATES  THAT  1!  HAS  BEEN  INITIALIZED 
BY  THE  USER  TO  A POSITIVE  DEFINITE  MATRIX. 

- VECTOR  OF  LENGTH  N CONTAINING  PARAMETER 

VALUES. 

ON  INPUT,  X MUST  CONTAIN  THE  INITIAL 
PARAMETER  r ST 1 MATES. 

ON  OUTPUT,  X CONTAINS  THE  FINAL  PARAMETER 
ESTIMATES  AS  DETERMINED  BY  ZXMIN. 

- VECTOR  OE  LENGTH  N*  (N+l)/2  CONTAINING  AN 

ESTIMATE  OF  THE  HESSIAN  MATRIX 
D**2F/(I)X(I)DX(.T) ) , I , J=1 , . . . ,N. 

H IS  STORED  IN  SYMMETRIC  STORAGE  MODE. 

ON  INPUT,  IP  IOPT  = 0,  ZXMIN  INITIALIZES  11 
TO  THE  IDENTITY  MATRIX.  AN  INITIAL  SETTING 
OF  H BY  THE  USER  IS  INDICATED  BY  I0PT=1. 

II  MUST  BE  POSITIVE  DEFINITE.  IF  IT  IS  NOT, 

A TERMINAL  ERROR  OCCURS. 

ON  OUTPUT, II  CONTAINS  AN  ESTIMATE  OF  11IE 
HESSIAN  AT  THE  FINAL  PARAMETER  ESTIMATES 
(1 • H. , AT  X ( 1) . X (2),...X(N) 

- A VECTOR  OF  LENGTH  N CONTAINING  AN  ESTIMATE 

OF  THE  GRADIENT  DI/DXfl),  I=1,...,N  AT  THE 
FINAL  PARAMETER  ESTIMATES.  (OUTPUT) 

- A SCALAR  CONTAINING  THE  VALUE  OF  THE  FUNCTION 

AT  THE  FINAL  PARAMETER  ESTIMATES.  (OUTPUT) 

- A VECTOR  OF  LENGTH  3*N  USED  AS  WORKING  SPACE. 

ON  OUTPUT,  WORK  (I),  CONTAINS  FOR  1=1,  THE  NORM 
OF  THE  GRADIENT  (1 .1:. , SQH‘1  (G(1)**2*G(2)**2+. . .«C.(N)**2) 
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PRECISION 

REQD.  IMSL  ROUTINES 
LANGUAGE 


1=2,  THE  NUMBER  OP  FUNCTION  HVALTIONS 

PERFORMED. 

1=3,  AN  ESTIMATE  OF  THE:  NUMBER  OF 
SIGNIFICANT  digits  in  the  final 
PARAMETER  ESTIMATES. 

ERROR  PARAMETER. 

JER  =0  IMPLIES  THAT  CONVERGENCE  WAS 

ACHIEVED  AND  NO  ERRORS  OCCURRED. 

TERMINAL  ERROR 

IER  = 129  IMPLIES  TH/AT  THE  INITIAL  HESSIAN 
MATRIX  IS  NOT  POSITIVE  DEFINITE.  THIS 
CAN  OCCUR  ONLY  FOR  IOPT  = 1. 

IER  - 130  IMPLIES  THAT  THE  INTERATION  WAS 
TERMINATED  DUE  TO  ROUNDING  ERRORS 
BECOMING  DOMINANT.  THE  PARAMETER 

' ESTIMATES  HAVE  NOT  BEEN  DETERMINED  TO 
NS1G  DIGITS. 

IER  - 131  IMPLIES  THAT  THE  INTERATION  WAS 
TERMINATED  BECAUSE  MAXFN  WAS  EXCEEDED. 

SINGLE 

UERTST, ZXMJN 

FORTRAN 


CALL  2XMIN  (FUNCT,N,NSIG, MAXFN, IOPT,X,H,G,F,W, IER) 


Purpose 


To  find  the  minimum  of  a function  f(x)  of  N variables  x = (x^  .Xj, . . • ,Xy) . It  is  assumed 
that  .the  gradient  vector  g(x)  = * •*  ■|>^)  and  t^lc  Hessian  matrix  H= 

exist,  but  the  user  is  not  required  tov supply  formulas  for  their  evaluation.' 

Algorithm 

ZXMIN  is  based  on  the  Harwell  library  routine  VA10A.  It  uses  a quasi-Newton  method. 

See  reference:  Fletcher,  R. , "Fortran  subroutines  for  minimization  by  quasi-Newton 
methods",  Report  R7125  AERE,  Harwell,  England. 

Programming  Notes 

1.  Accuracy  in  evaluating  the  function  f is  critical  when  highly  accurate  parameter 
values,  x-»  ate  required  (i.c.,  when  NSIG  is  greater  than  3 for  single  precision).  In 
these  casis,  it  is  advisable  to  evaluate  f in  precision  greater  than  working  precision 
(c.g.,  for  single  precision  AXMIN,  double  precision  can  be  used). 


2.  In  reference  to  significant  digit  tests  for  the  x vector,  leading  zeroes  to  the  right 
of  the  decimal  point  are  counted.  For  example,  both  123.45  and  0.00123  have  five  signif- 
icant digits.  Scaling  the  x vector  may  be  required  if  several  of  the  parameters  x.  arc 
much  less  than  1.0  to  obtain  the  same  number  of  significant  digits  for  each  of  the  x^ . 

3.  The  user  can  place  a limit  on  the  amount  of  computer  time  used  by  ZXMIN  by  setting 
MAXFN  appropriately.  In  general,  MAXFN=500  is  satisfactory  with  NS1G=3. 
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4.  The  accuracy  of  the  final  Hessian  matrix,  II,  is  not.  always  satisfactory.  Por  example 
when  a minimum  is  located  in  a very  few  steps,  the  estimate  may  be  inaccurate.  II  is  cal- 
culated by  updating  the  initial  Hessian  as  steps  occur. 

Example 

2 2 2 

Suppose  that  f=100(x2"Xj)  + (1-Xj)  is  ;to  be  minimized  and  parameter  estimates  Xj  and 

x 2 are  to  have  3 digit  accuracy. 

Input: 


EXTERNAL 

DIMENSION 

N=2 

NSIG=3 
MAXFN=500 
IOPT=0 
X(l)  = -1.2 
X(2)  = l.  0 


FUNCT 

\(2),  11(3),  G(2),  W(6) 


Call  ZXMIN (FUNCT, N, NS 1G.MAXFN, IOPT,X,H,G,F, W, IER) 

where  subroutine  FUNCT  has  the  following  code: 

SUBROUTINE.  FUNCT  (N,X,F) 

DIMENSION  X(N) 

F=100. * (X(2)-X(1)*X(1))**2+(1.-X(1))**2 
RETURN 


Output: 

IER  » 0 
X(l)  = .9999 
X(2)  = .9999 
10'5 


W(l)  10  (norm  of  gradient) 

W(2)  200  (number  of  evaluations) 

W(3)  3 (estimate  of  significant  digits  in  x) 


*.*  r:.v 


.n*V,  •*TS’  > 

: -y 


OPERATING  INSTRUCTION.0,  FOR  MF.XP  (RDF  6) 


INTRODUCTION 


This  manual  is  a guide  to  using  various  computer  programs  that  have 
application  in  linear  multivariable  systems  studies.  The  programs  are  described 
in  detail.  The  accompanying  FORTRAN  IV  listings  are  available. 

The  computational  algorithms  and  computer  programs  were  coded  by  the 
author  over  the  years  1967  - 1975.  They  range  in  scope  from  simple  matrix 
multiplication  routines  to  complex  routines  for  solving  matrix  Riccati  equations. 
The  coding  is  highly  efficient  with  considerable  attention  given  to  computa- 
tional speed.  Thus,  all  matrices  are  treated  internally  as  vector  arrays  in 
order  to  reduce  the  inherent  overhead  cost  of  dual  indexing. 

The  programs  have  been  used  successfully  for  several  years.  Their  worth 
and  efficiency  have  been  amply  demonstrated  in  numerous  applications. 


Labeled  COMMON  Blocks 

The  various  subroutines  require  certain  information  tc  be  made  available 
through  named  COMMON  blocks.  In  addition,  dummy  workspace  is  to  be  allocated 
in  the  MAIN  program.  The  blocks  required  are: 

1.  C0MM0N/MAIN1/NDIM,  NDIM1,  C0M1 

This  block  is  defined  in  the  main  pi'Ogram  where 

NDIM  = FORTRAN  specified  dimension  of  square 

arrays.  This  is  the  number  as  specified 
in  the  DIMENSION  declaration.  At  present, 

NDIM  must  be  80  or  less. 

NDIM1  = NDIM+1 

C0M1  = A dummy  workspace  (scratch -pad)  array. 

C0M1  must,  he  dimensioned  in  the  main 
program  to  be  of  dimension  NDIM  X NDIM. 


o 


2.  COMMON /MAIN2/COM2 


Defined  in  the  main  program.  COM2  is  an  array  containing  at  least 
NDIM  X NDIM  locations. 


3.  COMMON/ I NOU/ KIN,  KOUT,  KPUNCH 


Defined  in  the  main  program  for  input-output  control 


logical  unit  number  input  device 
logical  unit  number  output  device 
logical  unit  number  for  punch 


KPUNCH 


The  various  subroutines  generally  do  not  output  any  results  through  COMMON 
blocks.  The  transferring  of  data  and  computational  results  are  accomplished  via 
the  calling  sequences. 


Typical  first  few  cards  in  the  main  program  are 


DIMENSION  C0M1  (20,20),  COM2  (20, 20) .... (additional  arrays) 

C0MM0N/MAIN1/NDIM,  NDIM1 ,COMl/MAIN2/COM2/INOU/KIN, KOUT, KPUNCH 

NDIM=20 

NDIM1=21 

KIN=5 

K0UT=6 

KPUNCH=7 


The  maximum  allowable  dimension,  NDIM  = 80 


SUBROUTINE  MEXP 


(.  ) 


o 


Purpose  - Compute  the  NxN  matrix  exponential 
AT 

EA  «=  e 

using  a Chebyshcv  polynomial  approximation. 

Calling  Sequence  - CALL  MEXP  (N,A,T,EA) 

Inputs; 

N *=  size  of  A matrix  (NxN) 

A *=  parameter  array 
‘T  «*  scalar 

Outputs: 

EA  matrix  exponential 

Other  Programs  Used;  IDENT,  XNORM,  MMUL,  EQUATE,  DIAG 
Labeled  COMMON  block:  /MAIN1/ 

Computational  Algorithm:  The  Chebyshev  approximation 

n 


2x  = I (2)  + 2 £ I (2).T  (x) 
n m m 


m=l 


Wi1 


S ex1' 


i=0 

tVi  2x 

is  used  to  obtain  an  n = 12  order  approximation  to  e where 

I (2)  *»  modified  Bessel  function  of  the  first  kind 
m 


E rr 


k-0 


k!  (m+k) 1 
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(i) 


(2) 


(3) 


1 


o 


o 


o 


and 


T (x)  = Chcbyshev  polynomial 

in 


2s<m 


(-D 

2 s^O  ra_S 


SE 


m-s 

. 8 


<2x) 


m-2s 


C4) 


The  error  in  tlie  approximation  (2)  is  given  by 

|e  (x)j  < 21  (2)  • lx  I < 1 

I n I — n n 11  — 

Thus,  for  n *=  12,  an  accuracy  of 

|C12(X)I  — 2'3  X 

is  achievable,  theoretically. 

AX 

To  obtain  the  matrix  e , an  M^l  is  found  such  that 


(5) 


II A ||  T 

Ji 


< 1 


Eq.  (2)  is  then  used  to  obtain 

2.  « 


n 


2M  AT/2M_1 

i“0 

JI-1 


AT 


n f 


(6) 


, (7) 


AT  2 

Finally  e is  found  by  EA  = X . The  norm  of  the  resulting  error  matrix 

rv“9 


(theoretically)  satisfies  |e||  < 2.3  x 10  . 
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ZXMIN 

Minimized 


PRINT  L’s  & C’s 


START 


MAIN 


SUBROUTINE  ZXMIN 


SUBROUTINE  FUNCT 
CALL  MEXP 


Figure  16.  "Mini"  Flow  Chart  of  the  Parameter 
Estimation  Program 


,STCSi3.T7  8 0<*99,VANNOY,538  35 

t n O A O V C M — A C ft  - 


HON,T3  90,I07Q,CM6C  300  ,STCSt3.T7  80499,  VANNOY, 
FTN,R=3. 

ATTACH,  I MSI.,  IMSL, IO=LIDRARY  ,SN  = AS3. 

ATTACH  ,PFN»  CONTROL ,CY  = 4, IO=L720  0 33, SN-AFML. 
LIBRARY, IHSL,PFN. 

MAP, PART. 

0 


These  job  cards 
are  needed  for 
the  parameter 
estimator . 


PROGRAM  MIN 


Parameter  Estimation  Source  Listing 


74/74  OPT=l 


FTN  4. 6*446 


PROGRAM  MIN  (INPUT, OUTPUT) 
EXTERNAL  FUNCT 


THIS  PROGRAM  USES  A LEAST  SQUARES  FIT  OBJECTIVE  FUNCTION  MINIMIZING 
SUBRIOTINE,  ZXMIN,  TO  ESTIMATE  THE  INOUCTCR  A M J CAPICATOR  VALUES 
NEEDED  FOR  THE  E-LINE  TO  PRODUCE  A RECTANGULAR  PULSE  AT  THE  OUTPUT. 
THE  MAIN  PROGRAM  CONTAINS  THE  INITIAL  GUESS  FOR  THE  INDUCTOR  AND 
CAPACITOR  VALUES.  THE  MAIN  PROGRAM  ALSO  CONTAINS  PRINT  STATEMENTS 
AMO  VALUES  FOR  THE  VARIABLES  USED  BY  ZXMIN. 


DIMENSION  X(10) ,H(55) ,G(10) ,W(30) ,C(5) 
REAL  L (5 ) 


THE  FOLLOWING  VARIABLES  ARE  NEEOEO  BY  ZXMIN  (SEE  THE  OPERATING 
INSTRUCTIONS  FOR  ZXMIN). 


N * 10 
NSIG  = 12 
MAXFN  = 500 
IOPT  = 0 


THESE  VARIABLES  ARE  NEEOEO  BY  THE  MAIN  PROGRAM. 


THE  INITIAL  ESTIMATE  OF  THE  INDUCTOR  ANO  CAPACITOR  VALUES  WHICH 
GIVE  A RECTANGULAR  OUTPUT  VOLTAGE  PULES  ARE  GIVEN  BELOW! 

(THE  VALUES  ARE  IN  HENRYS  ANO  FARAOS  RESPECTIVELY) 


2.577E-07 
2.086E-07 
2. 17 lt-07 
2, 554F-07 
3.607E-07 
7.830L-06 


program  him 


74/74  OPT  = 1 


FTN  4.6+446 


,■  U 


C(2)  = 7.7826-06 
C(3)  = 8.521E-06 
C<4)  = 1. 079E-05 
C ( 5 ) = 2. 029E-05 


since  THE  parameter  values  used  dy  the  subroutine  zxmin  must  3E 
PUT  INTO  THE  X VECTOR,  THE  INOUCTOP  VALUES  BECOME  THE  OOO 
ELEMENTS  IN  THE  VECTOR  ANO  THE  CAPACITOR  VALUES  OECOME  THE  EVEN 
ELEMENTS  IN  THE  VECTOR. 


DO  1 1=1, NN 
NL  = 2*1-1 
NM  = 2*1 
X ( NL)  = L ( I) 
X(NM)  = C<I) 
1 CONTINUE 


CALLING  THE  SUBROUTINE  ZXHIN 


CALL  ZXHIN ( FUNCT , N,NSIG,MAXFN,IOPT,X,H,G,F,W,IER) 


PRINT  STATEMENTS 


PRINT  5 
5 FORMAT ( 1H1) 
PRINT*," 
PRINT*," 
PRINT*," 


PRINT*,' 

PRINT*,' 

PRINT*,' 

PRINT*,' 

PRINT*,' 


THE  INDUCTOR  ANO  CAPACITOR  VALUES  AS  CALCULATED  3Y 
ZXMIN  WHICH  GIVE, AS  CLOSE  AS  POSSIBLE,  THE  DESIRED' 
RECTANGULAR  PULSE  ARE)” 


00  10  1=1, NN 
IL  = 2*1-1 
IH  = 2*1 

PRINT*,”  L I ," ) = 
PRINT*,"  CC~,I,")  = 
10  CONTINUE 

PRINT*," 

PRINT*," 

PRINT*, "IER  = ", IER 
PRINT*,” 

PRINT*," 

PRINT*, "F  = ~,F 
PRINT*," 

PRINT*," 

00  15  1=1,3 


a,X(IL)  ,' 
•,X(IH)  ,’ 


G C", IL , ")  * ~,G(IL> 
G (“, IM, ")  = ", G ( IM) 


PROGRAH  MIN 
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PRINT*,"  W(",I,”)  = ” , W ( I ) 
15  CONTINUE 
STOP 
ENO 


- - ’'■•r 


c 


5 

C 

C 

C THIS  SUBROUTINE  TAKES  THE  INDUCTOR  AND  CAPACITOR  VALUES  OETERMINEO 
C BY  THE  SU  ROUTINE  ZX’IIN  AND  CALCULATES  THE  LEAST  SQUARES  OBJECTIVE 
10  C FUNCTION.  THE  FUNCTION  IS  FOUND  USING  NUMERICAL  INTEGRATION. 

C THROUGHOUT  THE  SUBROUTINE  LIMITS  HAVE  BEEN  SET  ON  CERTAIN  VARIABLES 
C SO  THAT  OVERFLOH  WILL  NOT  OCCUR. 

C 

C THE  COMMON  STATEMENT  AND  THE  VARIABLES  NOIM,  NOIM1,  ANO  COM1  ARE 
15  C NEEDED  3Y  THE  SUBROUTINE  WHICH  CALCULATES  THE  STATE  TRANSITION 

C MATRIX  CMEXP). 

C 

C 


20 


25 


30 


35 


40 


45 


50 


55 


C 

C 

C 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  FUNCT(N,X,F) 

DIMENSION  X(M> , A (11 , 11 ) , E A ( 1 i , 11) ,PAI (11,11) ,EAI (11,11) ,R(5) 
DIMENSION  COMKU.ll) 

REAL  J),J1,J2,J3,JT0T,K1 
COMMON/ MAIN1/NDIM,NOIM1,COM1 
NOIM  = 11 
NOIM1  = 12 


THE  FOLLOWING  IS  A LIST  OF  THE  VARIABLES  NEEOEO  IN  THE  LEAST 
SQUARES  FIT  OBJECTIVE  FUNCTION  SUBROUTINE. 

HI THE  DIMENSION  OF  THE  A MATRIX 


OELT THE  TIME  SU3INTERVAL  USEO  IN  THE  NUMERICAL  INTEGRATION 


THE  INITIAL  VOLTAGE  ON  THE  SOURCE  CAPACITOR 


TIME  WHEN  THE  OESIREO  PULSE  GOES  TO  ZERO  (SEE  EQUASION 
35  IN  THIS  TEXT). 


PULHT THE  OESIREO  PULSE  HEIGHT 


I PULSE- -TIME  WHEN  THE  OESIREO  PULSE  GOES  TO  ZERO, (TP),  OIVIOEO 
BY  BELT.  THIS  INTEGER  IS  AN  UPPER  LIMIT  ON  ONE  OF  THE 
SUMMERS  IN  THE  OBJECTIVE  FUNCTION  (SEE  EQASION  35  IN 
THIS  THESIS). 


TIME  WHEN  THE  OSCILLATIONS  OF  THE  ACTUAL  PULSE  SHOULO 
STOP,  -OTP),  OIVIOEO  BY  OELT.  THIS  INTEGER  IS  THE  UPPER 
LIMIT  ON  THE  SUMMER  IN  THE  OBJECTIVE  FUNCTION. 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


THE  FUNCTION  SUBROUTINE  USES  SEVERAL  OTHER  SUBROUTINES  THEY  ARE  I 

MEXP A SUBROUTINE  THAT  CALCULATES  THE  STATE  TRANSITION 

HATRIX. 

MMPY a SUBROUTINE  THAT  COMPUTES  THE  MATRIX  PROOUCT  C=A*B. 


C 

C 

C 

c 

c 

c 


NZERO— A SUBROUTINE  THAT  ZEROS  A MATRIX,  A=0.0. 

MIDENT-- A SUBROUTINE  THAT  CREATES  THE  IOENTITy  MATRIX. 

MEQUAL—  A SUBROUTINE  THAT  MAKES  ONE  SUBROUTINE  EQUAL  ANOTHER,  B=A 


74/74  0PT=1 


SUBROUTINE  FUMCT 
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JO  - 0,0 
J1  = 0.0 
J2  = 0.0 
J3  = 0.0 
K1  = 0.25 
F = 0.0 

CS  = 1.333E-04 
R( 1 ) = 15.0E-03 
RC  2)  = 15.0E-03 
R(3)  = 15.0E-03 
R(4)  = 15.0E-03 
R(5)  = 15.0E-03 
RL  = 0.165 


ALL  NON-ZERO  ELEMENTS  IF  THE  A MATRIX  ARE  OEFINEO  HER 


A(l,2)  = -1.0/CS 
A(2,l)  = 1.0/X(1) 

A ( 2| 2 > = -R  d>  /X  d) 

A 1 2 » 3 1 = -i.C/Xd) 

A(3,2)  = 1. 0/X (2 ) 

A ( 3 ,4)  = -1. 0/X (2) 

AU,3)  = 1.  0/X  (3) 

A ( 4 1 4)  = -R ( 2) /X l 3 ) 

A ( 4 ,5)  = -i.0/X<3> 

A ( 5 ,4 ) = 1.  0/X (4) 

A(5,6)  = -1 , 0/X (4) 

A ( b , 3 ) = 1.0/X15) 

A<6,6)  = -R ( 3) /X (5) 

A(6,7)  = -1.0/XT5) 

A<7,6)  = 1. 0/X (6) 

A<7,3>  = -1 , 0/X (6) 

A ( i ,7)  * 1.0/XC7) 

A(3,3>  - -R(4)/X<7) 

A<3 ,9)  = -1, C/X ( 7) 

At  3,3)  = 1, 0/X (8  > 

A ( 9 * 1 C ) = -1 . 0/X (8) 

A( ID) 9)  = 1 . C/X { 9) 

A ( 1 C > 1 0 ) = -R (5) / X (9) 
4(10,11)  = -1. 0/X (9) 

A(ll, 10)  = l.J/X(10) 

A(ll, 11)  = -1.0/ (RL*X(10) ) 


LIMITS  ARE  SET  HERE  ON  HOW  LARGE  THE  ELEMENTS  IN  THE  A MATRIX 
CAN  BE. 


EAMAX  = 610 . 0/DELT 
00  10  I« 1 , Ml 
00  10  J= 1 » N1 

IF(A(I, J) .GT. EAMAX)  GO  TO  11 
GO  TO  10 

11  A ( I , J)  a EAMAX 
10  CONTINUE 


SUBROUTINE  FUNCT 
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o 

175 

ISO 


165 


190 


195 


200 


205 


210 


215 


220 


225 


o 


c 

c«*m***«mm**m*mm«»<«**m**mm**mm**4«m«****m*m»m«m»*mm 

C************  * ********** 

c 

c . . 

c 

C THIS  SECTION  FINOS  THE  FUNCTION  F FOR  THE  ESTIMATED  VALUES 

C OF  THE  L‘S  AND  C*S.  NUMERICAL  INTEGRATION  IS  USED. 

C 

c 

c************. *********** 

c 

c 

Q # * * ♦ * ******  ******  *****  * 

C THE  SUBROUTINE  MEXP  IS  CALLED  HERE,  IT  CALCULATES  THE  STATE 

C TRANSITION  MATRIX  AT  TIME  T =0ELT . LIMITS  ARE  ALSO  SET  ON  HOW 

C LARGE  THE  MATRIX  CAN  BE.  THE  VARIABLES  USED  BY  MEXP  ARE  OEFINEO 

C IN  THE  OPERATING  INSTRUCTIONS  FOR  MEXP. 

Q*********************** 

c 

CALL  MEXP(N1,A,0ELT,EA) 

c 

OO  5 1=1, N1 
00  5 J= 1 , N1 

IF(EA<I,J).GT.1.0E*53)  GO  TO  6 
IF(EA(I,J) ,LT.-i.0E*53>  GO  TO  7 
GO  TO  5 

6 EA ( I, J)  = 1.0E+50 
GO  TO  5 

7 EA ( I , J ) = -l.GE*50 
5 CONTINUE 

C 

c ************  **********  * 

C r IS  THE  SOLUTION  TO  THE  OUTPUT  EQUASION  AT  THE  PARITCULAR 

C TIME  T=I*D£LT. 

c limits  are  also  set  here  on  hoh  large  y can  be. 

c*********************** 

c 

CALL  MIOENT (PAI,N1,N1) 

c 

00  100  I=6,IT0T 
CALL  MZER0(EAI,N1,N1) 

c 

CALL  MMPY(£A,PAI,cAI,Nl, N1,N1) 

C 

CALL  MEQUAL(PAI,EAI,N1,N1> 

c 

Y = EAI (Nl, 1)*X0 
C 

IF ( Y.GT • 1 • 0E*5C)  GO  TO  12 
IF(Y.LT.-1.0l*50)  GO  TO  13 
GO  TO  14 

12  Y = -1.0E+50 

PRINT*, “Y  IS  LESS  THAN  -1.0E*50“ 

GO  TO  14 

13  Y a 1.0E+5C 


96 


t 

c 

0 

e 

t 

b 

c 

c 

c 

c 

b 

c 
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2 30 


235 


240 


245 


250 


14 


15 


25 


20 


30 


PRINT* , "Y  IS  GREATER  THAN  1.0E+53' 
CONTINUE 

IFd.ER.l)  GO  TO  15 
IF(I.EQ.ITOT>  GO  TO  23 
IFU.GT.  IPULSE)  GO  TO  25 
PUL J = ( <PUIHT-Y)**2)*DELT 
J1  = J1  + PUL J 
PUL J =0.0 
GO  TO  30 
CONTINUE 

JO  = ( CPULHT-Y)**2)*OELT/2.0 

GO  TO  30 

CONTINUE 

ZEQOJ  = K1*(Y**2)*DSLT 
J2  = J2  + ZERO J 
ZEROJ  = 0.0 
GO  TO  30 
CONTINUE 

J3  = Kl*(Y**2)*DELT/2.0 
CONTINUE 


100 


Y = 0.0 

CONTINUE 

F = JO  + 

PRINT*," 

RETURN 

END 


J1  ♦ 
F = ' 


J2 

’,F 


+ J3 


1 


subroutine  mmpy 
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SUBROUTINE  MMPY(A,D,C,.1,  K,N) 


*****************  + ■***** 


THIS  SUBROUTINE  MILL  COB PUT  E THE  MATRIX  PROOUCT  C = A+B 


USAGE  1 

DIMENSION  A ( M , K ) , 3 ( K , N ) , C ( M , N ) 
CALL  MMPY(A,3,C,M,K,N) 


************** 


DIMENSION  A (M,K>  , 0 (K ,N ) , C (H , N ) 
DO  10  J=1,N 
00  10  1=1, H 
C(I,J)  = 0.0 
DO  10  L = 1 , K 

C(I,J)  = C( I,J>  * A<I,L) *B(L, J) 
10  CONTINUE 


00  20  1=1, H 

DO  20  J=  1 , N 

IP(C(I,J) .GT.1.0E+50)  GO  TO  35 
IF  (C(I,J) ,LT .-1, 0E+50)  GO  TO  30 
GO  TO  23 

30  C(I,J>  = -1.0E+50 
GO  TO  20 

35  C<I,J>  = 1.0E+50 
20  CONTINUE 


RETURN 

ENO 


******** 


SUBROUTINE  MZERO 
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a — 

e , 


c 

c 

c 

|c 

c 

r 

c 

c 

€ 

* 

c 

D 

c 

D 

c 

p 

p 

c 


o 


10 


13 


20 


25 


SUBROUTINE  MZERO ( A ,M,N> 

C 

c 

c*.  ***********  *i 

c 

c 

C THIS  SUBROUTINE  ZEROIS  A MATRIX  A=0.0 
C 

C USAGE  I 

C DIMENSION  A (M.N) 

C CALL  MZERO (A,M,N) 

C 

C 

c*«*«  ********* 


C 

c 


DIMENSION  A (M.N) 

00  10  1=1, M 
00  10  J=1,N 
A ( I , J ) = 0.0 
10  CONTINUE 
RETURN 
ENO 


99 
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SUBROUTINE  MIOENT 


SU3R0UTINE  MI DENT (D,M,N) 


C THIS  SU3R0UTINE  CREATES  THE  IDENTITY  MATRIX 


USAGE  * 

DIMENSION  D(M,N) 
CALL  MIOENT ( 0 , M, N) 


DIMENSION  D(M,N> 


CALL  MZER0(0,M,N) 


00  10  1=1, M 
0(1,1)  = 1.0 
10  CONTINUE 
GO  TO  30 


20  PRINT*, "ERROR*  MATRIX  IS  NOT  SQUARE  * CAN  *T  FORM  THE  IOENTITY  MAT 


30  CONTINUE 
RETURN 
END 


SUBROUTINE  MEQUAL 
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1 


5 


10 


15 


20 


25 


SUBROUTINE  MEQUAL  (A,B,M,N> 


C THIS  SUBROUTINE  MAKES  ONE  MATRIX  EQUAL  TO  ANOTHER  MATRIX 
C 

C USAGE! 

C DIMENSION  A tM  »N) ,9(H,H> 

C CALL  MEQUAL (A, n,M,N> 

C 

C 


C * * * 


C 

c 

01  MENS  I ON  A (Mi N) ,B(M,N) 
C 

00  10  1=1, M 
00  10  J = 1,N 
A ( I , J)  = 0(1, J) 

10  CONTINUE 
RETURN 
ENQ 


(J 


1 ) 


Typical  Output  from  the  Estimation  Program 


THE  IMJUCTOR  AND  CAPACITOR  VALUES  AS  CALCULATE!)  BY 
7XMIN  WHICH  GIVE, AS  CLOSE  AS  POSSIBLE,  THE  OESIREO 
RECTANGULAR  PULSc  ARE  I 


L«l) 
Ctl> 
L (2) 
C (2) 
L (3) 
C(3) 
L (4) 
C (4) 
L (5) 
C (5) 


.150479715151E-7 

POO  007  32&Clj8  30  0269 

.244114161551E-8 

0930C’779P3BS13154 

.51619J3590C1E-7 

OOOOG J5149 37-,99577 

.63937497369E-9 

OOCClC79ie663l904 

.270C40401654E-7 

00C3 20 23 671250329 


G ( 1)  «, 
G (2) 
GUI  * 
G (4) 

G (5)  * 
G(6) 
G (7)  = • 
G ( 3 ) = 
G (9)  = 
G ( 10) 


-727714059.3737 
a 43130141.9701 
50554195.6278 7 
= 35561723.86354 
997279413.9773 
= 2963141d. 78393 
•4.C97059552653E+102 
305234C5. 95086 
1.644976441214E+9 
a 16692738.39131 


IER  = 130 


F a 7655.104566558 


W(l>  = 4.097059552653E+102 
WC2)  = 418. 

W{3)  a -103.1792786124 


Vita 


Larry  W.  Vannoy  was  bom  on  20  December  195^  in  Black- 
well,  Oklahoma.  He  graduated  from  Hillcrest  High  School  in 
Dalzell,  South  Carolina  in  1973*  He  then  attended  the  Uni- 
versity of  South  Carolina  where  he  received,  in  1977 » the 
Bachelor  of  Science  Degree  with  concentration  in  Electrical 
Engineering.  Upon  graduation  he  received  a commission  in 
the  United  States  Air  Force.  In  September  he  was  directed 
to  the  School  of  Engineering,  Air  Force  Institute  of  Tech- 
nology, as  his  initial  active  duty  assignment. 


Permanent  Address:  5^13  Plantation  Drive 


Sumter,  South  Carolina  29150 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  fHTn-n  Oeta^Enlered) 

! REPORT  DOCUMENTATION  PAGE 


1.  REPORT  NUMBER 


AFIT/GE/EE/79-2 


14.  TITLE  (end  Subtitle) 


1?.  GOVT  ACCESSION  NO. 


DIGITAL  SIMULATION  AND  PARAMETER 
ESTIMATION  TECHNIQUES  FOR  THE 
E-LINE  PULSE  FORMING  NETWORK 


7.  AUTHORS 

Larry  W.  Vannoy 
2Lt,  USAF 

8.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Air  Force  Institute  of  Technology 
(AFIT-EN) 

Wright-Patterson  AFB,  Ohio  45433 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Air  Force  Aero  Propulsion  Laboratory 
Wright-Patterson  AFB,  Ohio  45433 


READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

3.  RECIPIENT'S  CATALOG  NUMBER 


5.  TYPE  OF  REPORT  A PERIOD  COVERED 


MS  Thesis 


6.  PERFORMING  O^G.  REPORT  NUMBER 


*•  CONTRACT  OR  GRANT  NUM8EROJ 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  & WORK  UNIT  NUMBERS 


12.  REPORT  DATE 

March  1979 

13.  NUMBER  OF  PAGES 

111 


14.  MONITORING  agency  NAME  a ADDRESSflf  different  from  Controlling  Office;  IS.  SECURITY  CLASS,  (ol  thie  report) 

Unclassified 

IS*.  OECLASSI  FI  CATION /DOWN  GRADING- 
SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (ot  thle  Report) 


Approved  for  public  release?  distribution  unlimited 


I 17.  DISTRIBUTION  STATEMENT  fof  the  *6sfr*cl  entered  in  Block  30.  II  dlllerent  from  Report) 


18.  SUPPLEMENTARY  NOTES 


16  MAY  i! 


Approved  for  public  release?  IAW  AFR  190-17  ^ '' wmi 

JOSEPH  P.  HIPPS.  Major.^^F 
Director  of  Information 

19.  KEY  WORDS  (Continue  on  reveram  aide  it  necmaamry  mnd  identity  by  block  number) 

Pulse-Forming  Networks 
Network  Simulation 
v Digital  Computer  Simulation 
Parameter  Estimator 

\ Line  Pulser  ■ 

ABSTRACT  (Continue  on  rarrran  aide  II  nacmaamry  mnd  Identify  by  block  number) 

A line-type,  voltage-fed,  pulse-forming  network,  called  the 
^E-line.^is  studied.  The  energy  to  be  delivered  to  the  load 
is  stored  in  a source  capacitor  at  the  front  of  the  network. 
Storing  all  the  energy  in  the  source  capacitor  leads  to  a more 
reliable  and  lightweight  system  than  is  achievable  by  a con- 
ventional network  where  all  the  capacitors  store  energy.  The 
network  is  to  deliver  100-200  rectangular  pulses  per  second  to 


FORM 
I JAN  71 


COITION  OF  I NOV  68  It  OBSOLETE 


SECURITY  CLASSIFICATION  or  THIS  PAGE  (T*h«i  Date  tottered) 


f 


f 1L. 

( Jf 


inmiASPiFiED 

security  Classification  of  this  PAGCfiWun  Dmim  F.nt*n<i) 

\ 

Block  20  continued. 

the  load  where  each  pulse  is  20-30  microseconds  in  duration  and 
has  a pulse  height  of  30,000  volts. 

A computer  program  is  developed  which  models  the  discharge  of  the 
network  and  is  used  to  analyze  the  output  pulse  shape  for  dif- 
ferent inductor  and  capacitor  values  within  the  network.  Another 
computer  program  is  developed  which  estimates  the  inductor  and 
capacitor  values  needed  for  the  network  to  give  a close  approxi- 
mation to  the  desired  rectangular  pulse.  The  simulation  program 
models  the  discharge  of  the  actual  network  very  well.  The  esti- 
i mation  program  as  applied  in  this  thesis  does  a poor  job  of  pre- 
l dieting  the  inductor  and  capacitor  values  needed.  The  desired 
rectangular  pulse  is  not  achieved.  Recommendations  are  included 
which  may  improve  the  simulation  program,  estimation  program, 
and  the  actual  pulse  shape . 
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